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Section  1 


SUMMMARY 

This  final  report  describes  work  performed  in  the  period  October  1991  to  September  1994  for  the  U.S. 
Air  Force  of  Scientific  Research  under  Contract  No.  F49620-91 -C— 0072.  Dr.  Julian  Tishkoff  was  the 
AFOSR  Project  Monitor,  The  work  was  performed  and  the  report  has  been  prepared  by  the  Mechanical  Sys¬ 
tems  Laboratory  at  the  GE  Corporate  Research  and  Development  Center,  Schenectady,  New  York.  This 
section  summarizes  the  technical  problem  which  was  considered.  The  conclusions  are  summarized  in  Sec¬ 
tion  3. 


Aeropropulsion  gas-turbines  generally  employ  non-premixed  turbulent  combustion.  These  en¬ 
gines  have  aggressive  design  objectives  -  including  supercruise,  increasing  thrust:weight  ratio,  low  ob¬ 
servables  in  the  exhaust,  higher  altitude  performance,  endothermic  fuel  capability,  and  low  pollutant  emis¬ 
sions  -  which  in  turn  lead  to  new  challenges  for  the  combustion  system:  decreasing  size,  increasing 
combustor  temperature  rise,  and  low  emissions  of  pollutants,  smoke  and  other  visibles.  Chemical  kinetics 
and  their  interactions  with  turbulence  control  these  phenomena.  Although  the  hypothesis  of  fast  chemistry 
or  "mixed-is-burned”  is  useful  for  understanding  traditional  design  issues  such  as  equilibrium  exit  tem¬ 
perature  fields  (“pattern  factor”),  it  cannot  address  these  challenges.  The  goal  of  this  three-year  research 
program  was  a  quantitative  understanding  of  turbuence-chemistry  interactions  in  the  above  systems. 

Computational  modeling  is  playing  an  increasingly  large  role  because  of  the  high  cost  of  prototype 
development  and  testing,  Today,  products  are  planned  and  sales  efforts  are  launched  on  the  basis  of  pre¬ 
dicted  performance,  and  metal  is  not  actually  cut  until  customers  are  lined  up.  Models  must  therefore  ad¬ 
dress  all  the  issues  of  concern  to  customers,  as  outlined  above.  This  is  a  more  ambitious  set  of  parameters 
than  the  traditional  requirement  of  only  the  pattern  factor. 

Because  the  principal  task  of  the  combustor  is  to  deliver  hot  combustion  products  to  the  turbine, 
mixing  alone  or  chemistry  alone  are  not  sufficient  criteria  by  which  to  judge  performance.  Mixing  and  chem¬ 
istry  are  simultaneous.  For  example,  the  large-scale  unsteadiness  sometimes  seen  in  combustors  is  often 
taken  as  evidence  that  the  in-flame  mixing  is  so  slow  that  it  becomes  process-limiting;  however,  there 
are  other  factors  at  play.  The  combustor  cannot  deliver  an  unsteady  flow  to  the  turbines  or  else  they  will 
be  damaged;  hence  the  exit  temperature  and  velocity  fields  must  be  steady  and  mixed.  Also,  much  of  the 
important  minor  species  chemistry  is  too  slow  to  be  affected.  The  Damkohler  numbers  (ratios  of  mixing  to 
chemical  times  for  all  important  mixing  scales  and  all  important  reactions)  are  an  excellent  way  to  character¬ 
ize  the  effect  of  mixing  in  the  context  of  chemistry.  Our  program  has  shown  that  in  aeropropulsion  combus¬ 
tors  most  important  Damkohler  numbers  are  in  the  distributed  reaction  regime.  This  includes  disparate  phe¬ 
nomena  such  as  flame  stability  (e.g.,  good  reproduction  of  blowout  limits)  and  CO  or  NOx  emissions. 

Our  research  program  addressed  the  technical  issues  by  the  development  (and  detailed  validation) 
of  computational  models  that  fjt  into  a  familiar  framework.  Otherwise  the  design  and  development  commu¬ 
nities  would  not  be  able  to  use  the  results.  It  also  seems  clear  that  different  modeling  approaches  will  be 
needed  to  address  different  issues.  Unsteady  methods  such  as  Large  Eddy  Simulation  or  Discrete  Vortex 
Models  may  be  used  to  isolate  features  of  the  flowfield,  such  as  mixing  layers.  Direct  Numerical  Simulation 
will  play  a  role  in  physical  sub-model  development.  These  methods,  however,  will  not  directly  transfer  in¬ 
formation  to  the  majority  of  engineers.  The  present  pdf  method  with  a  selected  reduced  chemistry  scheme 
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-  the  latter  having  been  tested  versus  full  chemistry  in  turbulent  and  not  laminar  combustion  -  is  superior  to 
others  in  the  context  of  gas-turbine  combustion. 

Statement  of  Work 

Taskt:  Literature  survey  of  well -characterized  bluff-body  stabilized  turbulent  non-premixed  flames.  Ex¬ 
periments  with  laser-based  spectroscopic  data,  particularly  Raman/CARS  and  fluorescence,  on  critical 
quantities  will  be  favored.  Compilation  of  data.  Recommendations  for  future  experiments.  (Leading  candi¬ 
dates  were  the  University  of  Sydney/SANDIA  “Reverse  Flow  Reactor,”  the  GE  bluff-body  stabilized  burner 
used  in  the  current  AFOSR  program,  the  ALTEX/SANDIA  burner  and  the  Air  Force  Aeropropulsion  Laborato¬ 
ry  burner.)  Recent  literature  on  piloted  jet  flames  will  also  be  surveyed  for  data  in  the  regime  of  intense  turbu¬ 
lence. 

Sections  2.1  and  2.2  describe  the  work. 

Task  2:  Formulation  and  discretization  of  equations  for  selected  turbulence-chemistry  sub-models,  in  el¬ 
liptic  form  appropriate  for  steady-state  multi-dimensional  computation.  The  sub-models  include: 

1 .  The  laminar  flamelet  model  with  a  model  for  scalar  dissipation  rate  and  with  full  chemistry. 

2.  Assumed  shape  pdf  method  with  lower  moments  from  balance  equations  and  with  reduced 
chemistry  schemes. 

3.  Direct  Monte-Carlo  pdf  evolution  methods  for  the  scalar  pdf,  with  reduced  chemistry  schemes. 

4.  Hybrid  partial  equilibrium-flamelet  approach,  which  treats  slow  (e.g.,  recombination)  reactions 
in  a  distributed  zone  mode  and  and  fast  (e.g.,  chain -branching)  reactions  in  a  flamelet  mode. 

5.  Semi-empirical  approach  which  treats  the  fine  scales  as  a  perfectly  stirred  reactor,  with  mass 
flux  into  this  regime  based  on  dimensional  analysis. 

6.  Models  which  may  emerge  as  research  progresses  at  GE  and  elsewhere. 

The  basic  flow  algorithm  can  be  modified  in  each  case,  as  required  to  cater  to  the  particular  sub¬ 
model,  but  in  general  will  have  second-order  (in  sense  of  stress/flux  correlations)  closure  on  turbu¬ 
lence  and  second-order  accurate  (in  sense  of  Taylor  series  truncation  error)  discretization. 

Picking  up  on  the  theme  that  mixing  and  chemistry  are  simultaneous,  a  new  tool  for  the  simulation  of  the 
interaction  between  (simultaneous)  prescribed  turbulence  and  full  chemistry  schemes  has  been  devel¬ 
oped.  We  argue  that  it  is  superior  to  other  tools  such  as  the  strained  laminar  flame  or  the  PSR.  It  also  traces  a 
direct  lineage  from  pdf/CFD  models  as  shown  below.  Sections  2.1  and  2.4  describe  the  work. 

Task  3:  Assessment  of  each  of  the  selected  turbulence-chemistry  sub-models  against  each  of  the  se¬ 
lected  datasets.  Comparisons  of  the  following  types  of  quantities  as  available  in  each  case: 

1 .  Mean  velocity  and  turbulence  fields. 

2.  Mean  major  species  and  temperature. 

3.  Fluctuations  (e.g.,  varianoe)  in  major  species  and  temperature. 

4.  Mean  minor  species. 

5.  Fluctuations  (e.g.,  variance)  in  minor  species 

6.  Important  covariances. 

Sections  2.2  and  2.5-2. 8  describe  the  work. 
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Section  2 

TECHNICAL  DISCUSSION 

Despite  significant  progress,  it  is  not  clear  how  turbulent  combustion  should  be  viewed.  Conceptual 
questions  arise  such  as  whether  turbulent  flames  act  as  ensembles  of  strained  laminar  “flamelets”  or  as 
broader  “distributed”  zones  of  reacting  species.  Engineering  models  of  combustors  rely  largely  on  simple 
approaches  such  as  fast  chemistry  (mixed  is  burned)  and  assumed  shape  probability  density  functions 
(pdf’s),  and  have  been  successful  in  terms  of  gross  features  such  as  the  profiles  of  combustor  exit  tempera¬ 
ture.  Present-day  demands  on  combustion  equipment  and  thus  on  computational  models  are,  however, 
increasing  the  need  for  advanced  understanding  of  turbulence-chemistry  interactions.  For  example,  (i) 
flameout  and  relight  in  turbine  combustors  are  related  to  interactions  with  chain-branching  reactions,  (ii) 
hydrogen  burnout  in  supersonic  combustors  is  related  to  interactions  with  recombination  reactions,  and  (iii) 
requirements  of  low  emissions  of  NOx,  CO,  smoke,  and  other  observables  are  related  to  non-equilibria 
among  species  such  as  oxyhydrogen  radicals  and  CxHy.  These  turbulence-chemistry  interactions  span 
several  orders  of  magnitude  in  Damkohler  number,  depending  on  the  specific  reaction  and  the  specific  tur- 
buience  process  in  question. 

New  models  for  turbulence-chemistry  interactions  in  air-breathing  combustion  have  been  devel¬ 
oped  and  compared  with  turbulent  flame  data.  Accomplishments  in  the  reporting  period  include: 

(1)  Development  of  the  hybrid  Monte  Carlo  joint  PDF/CFD  model  for  a  bluff-bodv  stabilized  flame.  Detailed 
comparisons  have  been  made  with  Raman  data  on  means  and  rms  of  mixture  fraction,  major  species  and 
temperature  in  non -premixed  bluff-body  stabilized  CO/H2/N2  and  CH4  flames; 

(2)  Development  of  the  “PaSR”  or  partially  stirred  reactor.  The  PaSR  offers  an  alternative  to  traditional  evalu¬ 
ations  in  laminar  flames,  which  are  not  relevant  to  turbulent  combustion.  Reduced  CO/H2  schemes  and 
reduced  methane  schemes  are  compared  with  full  schemes  in  the  PaSR  model.  These  studies  span  several 
decades  of  turbulent  mixing  frequency.  A  three -variable  reduced  scheme  for  complex  hydrocarbon  fuels  is 
compared  in  the  PaSR  with  the  starting  scheme  over  two  decades  of  turbulence  frequency.  Direct  compari¬ 
sons  of  different  mixing  models  are  conducted  in  the  context  of  a  fuN  chemistry  scheme. 

The  following  sub-sections  highlight  the  approach  and  the  results. 

2.1  Experimental  Set-Up  and  Hybrid  Monte  Carlo  Joint  PDF/CFD  Modeling 

The  pressure -corrected  mean  Navier-Stokes/assumed-shape  pdf/k-e  turbulence  model  does 
not  account  rigorously  for  the  turbulence-chemistry  interactions  of  dominant  interest,  but  affords  signifi¬ 
cant  geometric  flexibility  and  rapid  convergence  for  pressure-dominated  internal  flow  [Correa  and  Shyy 
(1987)].  On  the  other  hand,  the  joint  (velocity-composition)  pdf  transport  model  includes  turbulence  and 
chemistry  with  single-point  closure  [Pope  (1990)].  The  pdf  model  has  been  widely  used  to  compute  turbu¬ 
lent  jet  flames,  both  in  the  composition-only  form  [Chen  and  Kollmann  (1 989)]  and  in  the  joint  velocity- 
composition  form  [Pope  and  Correa  (1 986)]  .More  recently,  the  joint  pdf  model  has  been  extended  to  “ellip¬ 
tic”  (recirculating)  flow  and  applied  to  the  computation  of  CO/H2  bluff-body  flames  [Correa  and  Pope 
(1992, 1994)],  marking  a  step  toward  the  recirculation-stabilized  flowfields  found  in  practical  burners.  Bluff 
bodies  also  eliminate  the  pilot-stabilization  necessary  in  jet  flames  at  high  Reynolds  numbers. 

There  are  other  approaches,  intermediate  in  complexity  and  cost,  between  the  assumed  shape  pdf 
model  (#1)  and  the  joint  velocity- composition  pdf  transport  model  (#4).  The  "conditional  moment  clo- 
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sure”  model  (#2)  solves  conventional  time-averaged  field  equations  for  the  means  of  reactive  scalars 
conditioned  on  the  mixture  fraction  [Smith  et  al.  (1 992)] ,  and  is  applicable  only  when  fluctuations  about  this 
conditional  mean  are  small.  The  scalar  pdf  model  (#3)  does  not  treat  the  velocity  part  of  the  pdf,  using 
instead  conventional  turbulence  modeling  to  supply  the  scalar  (and  momentum)  fluxes  [Chen  etal.  (1989)]. 
Arguments  of  computational  cost  are  usually  made  to  support  models  #1-#3;  however,  the  speedup 
achieved  in  parallelizing  particle  tracking  pdf  computations  (accomplished  in  the  present  reporting  period 
of  this  contract)  shows  that  the  pdf  computations  are  tractable  in  practical  design  codes  [Correa  and  Braat- 
en  (1993)]. 

Methane  is  of  particular  interest  in  this  context:  (1)  for  scientific  reasons  because  it  affords  strong 
finite-rate  chemistry  effects  and  well-studied  reduced  kinetic  schemes,  and  (2)  for  practical  reasons  be¬ 
cause  there  is  a  significant  worldwide  natural  gas  economy.  For  example,  about  50  GW  of  new  gas-fired 
powerplants  are  being  sold  annually.  The  development  and  qualification  of  predictive  tools  (with  the  re¬ 
quired  geometry,  chemistry,  and  turbulence  capabilities)  will  aid  this  industry.  We  will  not  report  in  detail  on 
the  CO/H2  flame,  because  the  methane  flame  was  more  challenging.  The  CO/H2  flame  has  been  reported  in 
the  archival  literature  [Correa  and  Pope  (1992)]. 

Of  the  many  prior  bluff-body  flame  studies,  the  most  closely  related  is  that  of  Mash  et  al.  (1992)  who 
also  used  fluorescence-corrected  Raman  spectroscopy  of  methane-air  flames.  The  present  work  uses  a 
jet  which  is  twice  as  big  and  stresses  the  jet-dominated  regime,  which  is  less  susceptible  to  large-scale 
shedding  off  the  bluff  body;  greater  degrees  of  extinction  are  found.  Hence,  the  new  contributions  are  Ra¬ 
man  temperature  and  major  species  measurements  in  the  jet- dominated  regime  of  a  bluff-body  stabilized 
non-premixed  methane-air  flame. 

The  bluff  body  burner  and  the  inflow  conditions  are  shown  in  Fig.  1 .  The  axisymmetric  bluff  body  has 
an  outer  diameter  of  38. 1  mm.  with  an  axial  jet  of  diameter  “d”  3.18  mm.  located  in  the  center.  It  is  well  known 
that  parts  of  this  flowfield  (notably,  the  annular  shear  layer  shed  off  the  bluff  body)  can  be  dominated  by 
unsteady  effects.  Care  was  taken  to  operate  in  a  velocity  (jet  and  coflow)  regime  where  the  flame  was 
steady.  The  cold  jet  Reynolds  number  is  1 2,000,  based  on  the  jet  diameter  and  jet  exit  velocity  of  62.5  m/s. 
The  coflow  air  velocity  is  1 8  m/s.  The  back  surface  of  the  bluff  body  is  coated  with  a  thermal  barrier  material 
to  reduce  heat  loss.  The  flame  is  stabilized  by  the  recirculation  zone  provided  by  the  bluff  body.  The  tunnel 
cross  section  is  15  cm  x  15  cm,  large  enough  not  to  interfere  with  the  flame.  (The  calculations  are  made  for  a 
circular  duct  of  the  same  cross-sectional  area).  Visual  observations  of  methane  flames  at  various  air  and 
fuel  jet  velocities  were  used  to  select  the  conditions  tor  the  Raman-scattering  measurements  reported  be¬ 
low.  This  flame  is  anchored  by  the  bluff  body  and  is  almost  extinguished  in  the  neck  region,  before  re— ignit¬ 
ing  further  downstream. 

The  Raman  system  is  based  on  a  flashlamp  pumped  dye  laser  which  provides  pulses  of  —  1 J  in  2-4 
us,  within  a  0.2  nm  bandpass  at  488  nm  at  a  repetition  rate  of  10  Hz,  The  light  scattered  at  right  angles  is 
collected  by  two  lenses,  separated  in  frequency  by  a  3/4  m  spectrometer  and  is  detected  by  eight  photomul¬ 
tiplier  tubes.  The  photomultiplier  tubes  detect  anti-Stokes  vibrational  Raman  scattering  from  N2,  Stokes 
vibrational  Raman  scattering  from  N2,  O2,  H2,  H2O,  CO  and  CO2,  and  Rayleigh  scattering.  The  temporal 
resolution  (2-4  us)  of  the  technique  is  limited  by  the  laser  pulse  length,  the  spatial  resolution  (0.2  x  0.2  x  0.6 
mm)  is  limited  by  the  spectrometer  entrance  slit  and  the  collection  optics,  and  the  data  acquisition  rate  is 
limited  by  the  laser  repetition  rate.  The  flame  luminescence  was  broad  — banded  throughout  the  visible  re- 
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gion,  and  was  reduced  by  a  polarization  filter  in  the  collection  optics.  The  polarization  vector  was  aligned  to 
pass  horizontally  polarized  Raman  -  and  Rayleigh-scattered  light. 


Inlet  air  flow 
18  m/s 


Bluff  body 
38.1  mm  dia. 


Methane  jet 
3.18  mm  dia. 
62.5  m/s 


top  wall 
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recirculation  zone 
(sooty  yellow  flame) 
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neck  region 
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jet-like  flame 
(blue) 
30<x/d 


bottom  wall 

Figure  1.  Non-premixed  bluff- body-stabilized  methane  flame;  “d”  is  the  jet  diameter. 


The  instantaneous  temperature  on  every  laser  shot  was  determined  in  three  independent  ways: 

(1)  the  Stokes-anti -Stokes  (SAS)  ratio  from  N2,  which  yields  the  temperature  directly  [Drake  etal.  (1982)]; 

(2)  an  iterative  scheme  in  which  an  initial  temperature  is  guessed,  based  on  which  the  mole  fraction  of  all 
major  species  are  calculated  using  their  measured  vibrational  intensities.  The  mole  fractions  are  then  cor¬ 
rected  using  high-temperature  correction  factors  to  accountfor  changes  in  thefraction  ofthe  Raman  band 
falling  in  the  exit  slits  provided  for  the  respective  photomultiplier  tubes.  The  iteration  process  is  repeated 
until  the  sum  of  the  mole  fractions  is  unity;  and  (3)  Rayleigh  scattering;  Raman  data  on  the  major  species 
were  used  to  obtain  the  Rayleigh  cross  section  of  the  mixture  and  thus  provide  temperature  in  an  iterative 
manner.  The  last  two  methods  agreed  best,  to  within  10K  on  mean  temperature  and  to  within  50K  on  a 
shot— to-shot  basis.  Method  (2),  based  on  the  sum  of  mole  fractions  of  major  species,  is  reported. 

Initial  measurements  with  the  Raman  system  showed  that  there  was  significant  laser-induced 
fluorescence  (LIF)  interference  throughout  the  flame,  as  reported  elsewhere  [Masri  and  Dibble  (1989), 
Dibble  et  al.  (1990)].  The  LIF  was  fairly  broadband  and  contaminated  all  Raman  channels,  but  was  negligi¬ 
ble  in  the  Rayleigh  channel.  An  additional  difficulty  was  the  crosstalk  between  CH4  and  other  Raman  chan¬ 
nels,  primarily  O2.  For  488  nm  excitation,  the  Raman  interference  in  other  major  species  was  insignificant.  To 
account  for  these  two  additional  sources  of  contamination  in  the  signals,  the  system  and  calibration  proce¬ 
dures  were  modified.  Additional  photomultiplier  tubes  were  installed  in  two  Raman-free  regions  of  the 
spectra  to  monitor  the  LIF  on  a  shot-to-shot  basis.  These  PMT’s,  termed  FI  and  F2,  were  located  at 
540nm  (between  O2  and  CO2)  and  at  590nm  (between  CH4  and  H2O) ,  respectively.  These  two  signals  were 
found  to  correlate  very  well  with  each  other  and  with  all  other  Raman  signals,  so  that  the  use  of  FI  was  found 
to  be  sufficient  to  allow  corrections  in  all  other  channels. 
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A  calibration  procedure  was  used  to  correct  for  the  LIF,  following  the  procedure  of  Dibble  et  al.  (1990). 
A  38-nnm-dia.  honeycomb  burner  was  built  to  provide  laminar  diffusion  flames  of  30%CH4/70%CO.  The 
flame  was  visibly  sooty  and  yellow  at  the  downstream  end,  and  contained  enough  soot  precursors  and  LIF. 
The  calibration  factors  were  obtained  by  iteration.  First  the  raw  data  were  used  to  calculate  the  temperature 
(from  the  sum  of  mole  fractions)  and  mole  fractions  of  major  species.  The  contamination  caused  by  fluores¬ 
cence  was  particularly  evident  in  the  fuel-rich  regions  ofthe  flame.  The  next  step  was  to  estimate  correction 
factors  for  each  of  the  major  species  based  on  the  raw  data,  and  to  subtract  a  term  equal  to  the  product  of 
the  correction  factor  for  species  “i”  and  the  value  of  the  fluorescence  signal  measured  in  photomultiplier 
tube  FI .  The  calculated  temperature  and  mole  fraction  profiles  were  compared  with  predicted  values.  The 
process  was  iterated  to  convergence.  Figure  2  shows  the  temperature  and  selected  species  so  obtained. 
The  temperature  agrees  fairly  well  with  the  predicted  laminar  flamelet  calculations  [Chen  et  al.  (1989)]  for  an 
assumed  stretch  of  5/s;  the  calculations  did  not  depend  strongly  on  this  assumed  value.  The  corrections  for 
N2. 02  and  CH4  are  substantial,  and  the  corrected  species  data  again  agreed  with  the  laminar  flame  calcula¬ 
tions,  which  are  omitted  for  clarity.  Calibration  was  repeated  before  and  after  each  set  of  measurements. 


MIXTURE  FRACTION 


Figure  2.  Raw  and  fluorescence-corrected  Raman  data.  Temperature.  “X”=  raw  data;  =  cor¬ 
rected  data;  solid  line  =  laminar  flame  calculation  from  Ref.  5.  Species,  raw  data:  “♦”=  N2,  O2 
and  “X”=  CH4;  corrected  data:  “♦"=  N2,  O2,  and  “■”=  CH4. 
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Radial  profiles  of  temperature  and  mole  fractions  of  major  species  were  measured  at  axial  locations 
of  x/d  =  5, 10,  20,  30,  50  and  70  and  along  the  centerline.  Since  not  all  these  data  can  be  discussed  here, 
attention  will  be  confined  to  the  region  of  maximum  turbulence,  10  ^  x/d  ^  20. 

To  extend  the  pdf  model  from  parabolic  to  elliptic  flows,  an  iterative  pdf/elliptic  computational  fluid 
dynamics  (CFD)  approach  has  been  developed.  The  highly  non-uniform  CFD  grid  accounts  for  the  dispa¬ 
rate  flow  scales  imposed  by  the  jet  and  bluff  body  dimensions,  with  a  second-order  accurate  numerical 
discretization  scheme  to  eliminate  numerical  “diffusion”  errors  [Correa  and  Shyy  (1987)]. 

The  transport  equation  for  the  joint  pdf  P0Ey;x,t),  which  describes  the  joint  probability  of  the  fluid 
composition  being  the  vector  \|2  and  the  fluid  velocity  being  the  vector  V  (at  the  position  x  and  time  t)  is: 

9(U-)  f  +  9«l-)  u, 


where  V  is  the  velocity  field,  the  density  q  is  a  function  of  the  scalar  variables  denoted  by  the  vector  and 
the  reaction  rate  is  S<i(  it).  The  six  terms  describe,  respectively,  the  evolution  of  the  pdf;  the  transport  of  the 
pdf  in  physical  space  (x;);  the  transport  of  the  pdf  in  velocity  space  Vj  by  body  forces  and  the  mean  pressure 
gradient;  the  transport  of  the  pdf  in  composition  space  tjia  by  reaction  Sa(  it);  the  transport  of  the  pdf  in 
velocity  space  Vj  by  viscous  stresses  and  the  fluctuating  pressure  gradient;  and  the  transport  of  the  pdf  in 
composition  space  by  molecular  fluxes.  A  one-point  statistical  description  in  terms  of  the  joint  pdf  of  the 
velocity  and  these  scalars  is  sought.  If  the  flow  is  statistically  stationary,  all  one- point  statistics  depend  only 
on  the  spatial  coordinates.  All  one-point  statistics  are  recovered  from  this  pdf  because  the  composition 
is  a  known  function  of  the  above  scalars.  The  velocity-composition  joint  pdf  evolution  model  relaxes  many 
of  the  assumptions  made  in  the  standard  closure,  such  as  pdf  shape,  statistical  independence  of  scalars, 
and  gradient  diffusion  of  scalars.  Closure  of  the  non-linear  chemical  source  term  or  “turbulent  fluxes”  of 
scalars  are  given  directly  by  the  pdf. 

Discussing  specifically  the  case  of  methane,  which  is  more  complex  than  CO/H2),  the  chemistry  is 
described  by  the  mixture  fraction  5,  and  the  reactive  scalars  in  the  system.  In  the  pdf  transport  equation, 
these  become  the  independent  variables  denoted  by  V,  i[i,  and  {|)k  (k=1,...4)  respectively.  The  joint  pdf 
evolves  in  this  eight-dimensional  velocity-composition  space  as  well  as  in  the  two-dimensional  x-r 
physical  space. 

Turbulent  mixing  of  the  scalars  is  modeled  as  a  linear,  deterministic  relaxation  to  the  local  mean, 
sometimes  called  Interaction-by-Exchange-with-the-Mean  (lEM)  [Borghi  (1988)].  The  lEM  model  has 
been  isolated  and  studied  in  detail  in  the  Partially  Stirred  Reactor  model,  described  next,  and  compared  with 
other  mixing  models.  Regarding  the  velocity  term,  the  fluctuating  component  of  acceleration  (arising  from 
the  fluctuating  pressure  gradient  and  viscous  forces)  is  modeled  by  the  simplified  Langevin  model  [Pope 
(1985),  Pope  (1990)].  The  Seshadri-Peters  (1990)  4-step  reduced  scheme  is  adopted  for  the  kinetics.  A 
look-up  table  of  reaction  rates,  density  and  temperature  is  constructed  on  a  non-uniform  20x10x10x10 
X  10  grid. 

The  axisymmetric  elliptic  mean  flow  CFD  model  and  the  pdf  model  communicate  with  each  other 
iteratively.  On  each  time  step,  the  fields  of  the  mean  velocity  and  the  turbulence  frequency,  obtained  from 
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the  local  turbulence  kinetic  energy  (tke)  and  dissipation  rate,  are  passed  from  the  CFD  model  to  the  pdf 
model.  A  shift  and  a  uniform  stretching  in  V  space  are  applied  to  the  pdf  so  that  the  mean  velocity  and  the  tke 
of  the  two  sub-models  are  in  agreement.  Thus  the  no-slip  boundary  condition  is  automatically  satisfied  at 
walls.  Stochastic  Lagrangian  particle  evolution,  per  the  lEM  model  and  the  reduced  scheme,  occurs  at  the 
above  frequency.  The  mean  density  field  is  passed  back  to  the  CFD  model,  and  the  two  sub-models  are 
iterated  to  convergence.  The  75  x  60  cell,  =  1 0^  particle  calculation  was  run  until  a  statistical  steady-state 
was  achieved.  For  display  purposes,  lower  moments  such  as  the  means  were  averaged  over  the  last  200 
time-steps  of  the  pdf  evolution,  reducing  statistical  fluctuations  in  the  results.  Scatterplots  were  prepared 
from  the  calculations  by  saving  several  realizations  of  the  pdf  particle  array  after  stochastic  convergence, 
and  then  post-processing  the  more  than  10®  particles  so  obtained. 

2.2  Comparison  with  Raman  Data  from  Non-Premixed  Bluff-Body  Stabilized  Flames 

The  following  discussion  focusses  on  the  region  of  strongest  turbulence  (10^x/d^20),  where  the 
off-axis  recirculation  zone  provides  intense  mixing  at  the  edges  of  the  jet. 

Computed  radial  profiles  of  the  mean  and  r.m.s.  mixture  fraction  at  x/d=20  compare  well  with  the 
Raman  data  (Fig.  3).  The  measured  radial  profiles  are  shown  in  full,  revealing  the  degree  of  symmetry,  while 
the  computed  profiles  are  by  assumption  axisymmetric.  The  comparisons  indicate  that  the  calculated  mix¬ 
ture  fraction  field  is  accurate  enough  to  permit  a  meaningful  evaluation  of  the  reactive  quantities.  The  mean 
temperature  peaks  in  the  recirculation  zone  with  a  maximum  of  less  than  1000K,  a  result  of  the  strong  turbu¬ 
lence.  Computations  and  data  agree  quite  well  (Fig.  4),  although  the  temperature  is  underpredicted  by  al¬ 
most  250K  along  the  centerline.  The  mean  major  species  profiles  also  agree  quite  well  (Fig.  5).  The  mean  O2 
is  depleted  in  the  wake  of  the  bluff— body  but  coexists  with  mean  CH4,  a  consequence  of  finite— rate  chemis¬ 
try.  The  O2  returns  to  ambient  levels  at  the  edge  of  the  bluff-body  (rs0.02  m).  The  model  predicts  more 
mean  O2  than  measured  at  the  centerline,  which  is  consistent  with  underprediction  of  the  mean  tempera¬ 
ture  in  Fig.  4. 

Scatterplots  provide  an  instructive  format  in  which  to  study  turbulence— chemistry  interactions,  and 
best  utilize  the  power  of  the  time-  and  (not  quite)  space -resolved  Raman  spectroscopy  and  the  pdf  model. 
The  measured  temperature-mixture  fraction  scatterplot  using  all  data  from  x/d=10  and  x/d=20  is  shown 
in  Fig.  6  (a),  along  with  the  calculated  laminar  flame  profile  for  a  stretch  of  5/s  [Chen  et  al.  (1988)].  Unlike 
in  the  CO/H2  flame  studied  previously,  bimodality  is  clearly  evident.  A  significant  number  of  points  has  mix¬ 
ture  fraction  values  close  to  stoichiometric  but  temperatures  which  manifest  localized  extinction,  while  other 
points  are  clustered  along  the  line  of  strained  flamelet  temperatures.  The  physical  picture  that  emerges 
from  the  scatterplots  is  that  of  a  recirculation  zone  (which  anchors  the  turbulent  flame),  connected  to  a  jet— 
like  region  by  a  narrow  neck  of  high  shear  with  a  significant  amount  of  local  extinction. 

The  corresponding  calculated  scatterplot  (Fig.  6  (b))  is  similar,  although  the  predicted  points  exhibit 
a  greater  trend  towards  extinction.  On  a  related  note,  stochastic  calculations  of  lean  methane— air  combus¬ 
tion  in  a  Partially  Stirred  Reactor  [12]  showed  that  the  parent  starting  scheme  used  to  develop  the  4— step 
model  agreed  with  a  77-step  scheme  reasonably  well  when  the  mean  temperatures  were  above  1500K, 
but  prematurely  predicted  blowout  relative  to  the  77— step  scheme.  This  behavior  is  similar  to  the  overpre¬ 
diction  toward  extinction  in  the  present  study. 
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Figure  5.  Comparison  of  mean  CH4  and  O2  profiles  at  x/d=20= 

The  coexistence  of  fuel  and  oxygen  is  again  apparent  in  Ych4-Yo2  scatterplots  (not  shown)  from 
both  the  data  and  the  model.  The  measured  and  calculated  scatterpoints  at  x/d=20,  both  of  T-^  and 
Ych4“Yo2.  were  closer  to  the  chemically  “frozen”  line  than  those  at  x/d  =  10,  in  accordance  with  Raman 
data  and  with  visual  observations  of  the  flame. 

The  calculated  CO  scatterplots  have  maxima  of  3%  (inset:  Figure  7)  whereas  the  data  peak  at  about 
10%  (Figure  7),  well  above  the  flamelet  maxima.  Similar  10%  levels  were  measured  in  Masri  et  al.  (1992) 
bluff  body  flame  and  in  pilot-stabilized  flames  [Chen  et  al.  (1989)],  and  2-3%  peaks  were  predicted  in  the 
latter  using  the  4-step  scheme  within  a  (scalar)  pdf/Reynolds  stress  model.  Hence,  this  discrepancy  on 
CO  maxima  has  appeared  in  diverse  circumstances  (but  always  In  combustion  gases  that  are  near  local 
extinction).  There  are  many  potential  contributors,  including:  (i)  the  assumption  of  a  steady  state  for  the 
radicals  in  the  4-step  mechanism;  (ii)  the  errors  in  Raman-based  CO  and  CO2  data,  as  discussed  above, 
and  as  seen  by  direct  comparison  with  predicted  mean  CO  and  CO2  profiles  in  the  CO/H2  bluff  body  flame 
of  Correa  and  Pope  (1992);  and  (iii)  neglect  of  phenomena  such  as  unsteady  flamelets  or  micromixed  gases 
(perfectly  stirred  reactors),  which  have  been  shown  to  lead  to  high  CO  [Maup  et  al.  (1990),  Chen  and  Dibble 
(1991)].  In  intense  turbulence,  however,  the  microscale  may  be  better  simulated  by  the  PaSR  -  described 
in  the  next  section  -  since  it  is  the  degenerate  form  of  the  pdf  equation  for  spatially  homogeneous  systems. 
For  example,  at  30  atm  the  PaSR  indicates  that  approximately  2%  peak  CO  levels  are  encountered  until  the 
fuel  is  pyrolyzed  and  CO  oxidation  can  commence  [Correa  (1 994)] .  It  seems  clear  that  CO  Is  an  important 
clue  to  the  microstructure  of  highly  turbulent  combusting  gases,  and  that  accurate  measurements  will  be 
critical. 
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MEASURED  MASS  FRACTION  OF  CO,  Y{CO) 


MEASURED  MASS  FRACTION  OF  CH4.  Y(CH4) 

Figure  7.  CO  -  methane  scatterplot  at  x/d  =  10  and  x/d=20. 

Because  the  joint  pdf  contains  the  velocity  components  of  each  particle,  the  scalar  fluxes  can  be  ex¬ 
amined.  Figure  8  shows  the  computed  radial  transport  ^  of  the  mixture  fraction  at  x/d=20,  calculated  by 
summing  over  particles  in  twenty  radial  bins.  Also  shown  is  the  gradient  diffusion  flux  -(m-i/q)  the  latter 
computed  from  the  mean  fields;  is  the  turbulent  viscosity  computed  from  the  k~e  equations.  Both  fluxes 
were  calculated  a  posteriori,  since  they  were  not  needed  for  the  main  computation.  Figure  8  shows  that 
in  the  presence  of  the  strong  radial  gradient  3|/3r,  (radial)  transport  of  the  mixture  fraction  Is  consistent  with 
the  notion  of  gradient  diffusion,  i.e.,  there  is  no  counter— gradient  diffusion.  The  magnitude  of  v'^'  is,  howev¬ 
er,  much  greater  than  that  of  —(ut/g)  0f/3r.  The  indicated  "turbulent  Schmidt  number"  is  about  0.4. 
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TRANSPORT  OF  CONSERVED  SCAU\R,  m/s 


Figure  8.  Radial  flux  of  conserved  scalar  at  x/d=20. 

The  study  demonstrates  several  points: 

1 .  The  bluff  body  burner  provides  a  strongly  turbulent  field  leading  to  localized  extinction,  without  the  need 
for  a  pilot  flame.  Thus  the  two -stream  nature  of  the  problem  is  preserved,  unlike  many  piloted  jet  flame 
studies  where  the  composition  or  the  excess  enthalpy  of  the  pilot  flame  can  cause  modeling  difficulties.  The 
present  recirculation-stabilized  flame  is  also  much  closer  to  practical  burners. 

2.  By  correcting  for  fluorescence,  Raman  measurements  can  be  made  in  bluff-body  CH4-air  flames;  how¬ 
ever,  the  errors  in  certain  species  (e.g.,  CO  and  CO2)  may  be  so  large  that  models  should  not  be  changed 
on  the  basis  of  such  data  alone. 

3.  Given  the  similarity  in  scatterplots,  it  is  clear  that  the  pointwise  structure  of  the  above  bluff  body  flame 
and  the  piloted  jet  flame  are  quite  similar.  A  greater  degree  of  local  extinction  is  measured  here. 

4.  The  limitations  of  pdf  shape  assumption,  statistical  independence  of  scalars,  and  gradient  diffusion,  are 
removed  in  this  model.  The  consequences  can  be  seen  in  joint  scatterplots  and  in  convective  fluxes. 

The  acquisition  of  Raman  data  and  the  3-velocity/5-scalar  joint  pdf  calculation  in  this  bluff  body 
methane  flame  takes  each  "discipline”  to  the  limits  of  the  present  state  of  the  art.  Any  significant  further 
progress  will  require  improvements  in  major  species  measurements,  complementary  velocity  and  minor 
species  measurements,  and  parallel  computers.  Reduced  chemistry  schemes  that  relax  steady  state  as¬ 
sumptions  (and  are  likely  to  require  additional  scalars)  will  have  to  be  developed  and  assessed  in  simpler 
contexts. 
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2.3  The  Need  for  an  Improved  Prototype  of  Highly-Intense  Turbulent  Combustion 

The  laminar  flame  speed  depends  on  the  fuel-air  equivalence  ratio  “O"  and  the  pressure  “p,”  and 
has  been  computed  for  the  methane-air  system  by  several  research  groups  with  reasonable  agreement 
[Smooke  and  Giovangigli  (1990),  Seshadri  and  Peters  (1990)].  Taking  the  inlet  temperature  Tjn  to  be  300K 
in  the  base  case,  the  flame  speed  drops  from  about  40  cm/s  in  the  stoichiometric  atmospheric-pressure 
case  to  about  10  cm/s  in  lean  atmospheric-pressure  flames,  and  similarly  in  stoichiometric  high-pressure 
flames.  The  flame  speed  would  decrease  further  with  strain.  On  the  other  hand,  recent  calculations  by 
Smooke  (1993)  indicate  that  the  flame  speed  rises  with  inlet  temperature;  this  behavior  is  relevant  to  high- 
pressure  combustion  because  the  fuel-air  premixture  is  preheated  by  the  (usually  nearly  adiabatic)  com¬ 
pression  process.  It  is  clear,  however,  that  the  laminar  methane-air  flame  speed  is  less  than  2  m/s  under 
conditions  Tjp  <  700K  and  O  <  0.6.  The  latter  are  typical  of,  for  example,  a  low-emissions  lean  premixed 
gas-turbine  combustor  operating  at  a  pressure  of  15  standard  atmospheres  in  the  burner. 

Next,  consider  the  above  flame  speed  relative  to  the  turbulence  intensity.  Taking  a  typical  combustor 
mean  flow  velocity  of  100  m/s  and  a  typical  turbulence  intensity  of  10%.  the  turbulent  velocity  fluctuations 
are  on  the  order  of  10  m/s,  i.e.,  on  the  order  of  ten  times  the  laminar  flame  speed.  The  topology  of  the  laminar 
flame  will  not  survive  in  turbulence  of  such  intensity.  This  argument  can  also  be  advanced  by  comparing 
the  time-scales  of  all  relevant  mixing  processes  (given  by  the  turbulence  spectrum)  with  the  time-scales 
of  all  relevant  chemical  reactions.  The  ratios  are  called  the  Damkohler  numbers:  values  much  less  than  unity 
suggest  behavior  approaching  stirred  reactors  [Bilger  (1 988)] ,  while  values  much  larger  than  unity  are  nec¬ 
essary  for  the  flamelet  model  to  hold  [Peters  (1 986)] .  It  can  easily  be  shown  that  in  real  systems  the  Damkoh¬ 
ler  numbers  span  several  orders  of  magnitude  about  unity,  so  flamelet  behavior  does  not  exist  for  all  the 
reaction/mixing  process  pairs  that  are  operative  [Correa  (1992)]. 

From  the  above,  it  may  be  concluded  that  an  intense  combustion  process  should  not  be  viewed  as 
the  motion  of  a  high-speed  flamefront.  Another  viewpoint  is  that  of  chemical  reactions  in  a  highly-loaded 
reactor  devoid  of  spatial  features.  The  “linear  eddy”  model  provides  a  one-dimensional  framework  span¬ 
ning  these  two  extremes  [Kerstein  (1991)].  Transport  by  turbulent  “eddies”  is  accounted  for  by  random 
rearrangement  of  the  scalar  fleld(s)  along  a  line,  at  a  rate  and  with  (a  distribution  of)  length  scales  in  accor¬ 
dance  with  fundamental  scaling  rules  taken  from  the  theory  of  homogeneous  turbulence.  In  between  these 
folding  events,  molecular  diffusion  (and  chemical  reactions)  can  occur  along  the  line.  As  the  turbulence 
amplitude  and  frequency  increase,  the  relevance  of  reaction -diffusion  interfaces  (flamefronts)  decreases. 
Ultimately,  the  velocity  field  and  the  spatial  structure  become  less  Important  and  a  stirred  reactor  description 
would  be  more  appropriate. 

2.4  Development  of  the  Partially-Stirred  Reactor  Model 

A  stochastic  model  operating  in  Ns-dimensional  composition  space  has  been  developed  forthe  lat¬ 
ter  class  of  spatially  homogeneous,  but  not  mixed,  reacting  flows.  The  model,  called  the  Partially  Stirred 
Reactor  or  “PaSR,”  Is  described  next. 

The  PaSR  model  has  two  applications: 

(i)  as  the  closure  model  for  the  multi-dimensional  joint  pdf  model.  In  this  application,  simplified  chemical 
schemes  are  typically  used.  An  example  is  the  combined  CFD  (Computational  Fluid  Dynamics)  and  joint 
pdf  calculation  of  a  bluff-body  (recirculation)  stabilized  CH4  flame,  which  used  five  scalar  variables  to  de- 
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scribe  the  chemistry;  and 

(ii)  for  standalone  calculations  of  the  intensely  turbulent  but  not  perfectly  stirred  regime,  decribed  in  the  fol¬ 
lowing  sections. 

The  PaSR  was  developed  in  this  program  during  this  reporting  period  [Correa  (1993  a)].  The  govern¬ 
ing  equations  may  be  derived  from  the  joint  pdf  equation  [Pope  (1990)]  by  specialization  to  homogeneous 
but  unmixed  flow,  obtaining: 


iE  +  '>[<fu>pi 

St  SiJ). - (2) 

Physically  speaking,  reactants  flow  into  a  "box”  whose  contents  -  reactants,  products  and  inter¬ 
mediates  -  are  mixed  by  turbulence  of  a  prescribed  frequency  (Fig.  9).  A  steady  state  of  “unmixedness” 
is  maintained  by  the  compositional  difference  between  the  inlet  stream  and  the  contents  of  the  reactor.  The 
pdf  P(Yk)  of  mass  fractions  of  species  in  the  reactor  is  represented  by  the  Np-  particle  ensemble 

Y|;^',Yf> . Y<"' . Y[^P>:  k  =  1 . N, 

where  Ng  =  number  of  species  and  Np  =  number  of  particles.  Scalar  mixing  is  accounted  for  by  the  “Inter- 
action-by-Exchange-with-the-Mean"  (lEM)  submodel  [Borghi  (1988)].  The  particle  equations  are 


M 

=  -C,o.)(Y^Y,)  +  W<"'^:  k  =  1 . N3  (3a) 

where  w  is  the  mixing  frequency  in  the  lEM  model,  Wk<")  is  the  molar  production  rate  of  species“k”  per  unit 
volume  for  the  particle,  is  the  molecular  weight  of  species  "k”  and  is  the  density  of  the  n*^  particle. 
Equation  3  (a)  is  solved  as  shown,  without  fractional  steps;  therefore,  mixing  and  chemistry  are  accounted 
for  simultaneously.  The  corresponding  equation  for  particle  temperature  is 


dt 


dH(p> 

dt 


Ns 

Sh, 

n  =  1 


(3b) 


where  is  the  total  enthalpy  of  particle  “n.”  Hence,  the  PaSR  is  described  by  a  coupled  system  of  (Ng+I) 
X  Np  first-order  stiff  ordinary  differential  equations  (o.d.e.’s). 


The  mass  flow  rate  m  into  the  PaSR  is  discretized  into  Nin  particles  per  global  time  step  (Fig.  9).  Par¬ 
ticles,  selected  randomly  from  the  ensemble,  flow  out  at  the  same  mass  flow  rate.  Particles  in  the  ensemble 
evolve  according  to  Eqs.  3.  The  integration  is  continued  until  a  stochastic  steady  state  is  achieved.  Scatter- 
plots,  pdf’s,  and  correlations  of  interest  can  be  obtained  from  the  steady  state  ensemble.  The  global  resi¬ 
dence  time  can  be  computed  as  t  =  gV/m,  where  q  is  the  mean  density  (obtained  by  appropriately  sum¬ 
ming  over  particle  masses  and  volumes)  and  V  is  the  reactor  volume;  x  also  converges. 
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ENSEMBLE  AVERAGES 

Np  =  number  of  particles  in  ensemble 
p  =  mean  density 

=  mean  species  mass  fractions 
T  =  mean  temperature 


Figure  9.  The  PaSR  model.  The  PaSR  may  be  viewed  as  a  single  computational  cell  in  the 
multi-dimensional  CFD/pdf  model  and  hence  previews  the  eventual  capability  of  the  latter. 

2.5  Parallelization 

Expensive  numerical  integration  techniques,  such  as  the  Gear  method,  are  required  because  of  the 
stiffness  of  the  o.d.e.’s  when  complex  kinetics  are  considered.  Therefore  the  PaSR  model  entails  consider¬ 
able  computational  cost.  Further  examination  indicates  that  solving  the  systems  of  ordinary  differential 
equations  for  each  of  the  particles  (i.e.,  Eqs.  2  and  3)  takes  up  almost  all  the  computational  time.  As  a  pre¬ 
lude  to  the  present  study,  the  PaSR  algorithm  was  evaluated  on  parallel  computers  with  a  view  to  reduction 
of  run  time.  The  term  "parallel”  refers  to  a  computer  architecture  in  which  multiple  linked  processors  work  in 
parallel  on  different  particles,  according  to  techniques  described  below.  The  lEM  model  lends  itself  ex¬ 
tremely  well  to  parallel  computing  because  the  terms  in  the  (Ng  +  I)  x  Np  particle  equations  contain  proper¬ 
ties  of  only  (i)  the  present  particle  and  (ii)  the  ensemble;  it  is  not  necessary  for  any  one  particle  to  "know” 
about  other  particles  individually,  but  only  indirectly  through  the  interaction  with  the  mean  calculated  from 
the  previous  time  step.  As  a  result,  all  particles  can  be  advanced  in  parallel  within  a  time  step.  Individual 
particles  or  groups  of  particles  can  be  computed  on  different  processors.  Parallelization  of  this  sort  is  re¬ 
ferred  to  as  particle  partitioning,  since  each  processor  operates  on  a  sub-ensemble  of  particles.  At  the  end 
of  the  time  step,  communication  between  the  processors  is  required  to  allow  the  mean  properties  to  be 
computed  for  use  in  the  next  time  step. 

The  approach  taken  was  to  parallelize  the  PaSR  algorithm  in  as  simple  a  fashion  as  possible.  Since 
the  solution  of  the  particle  equations  takes  up  almost  all  of  the  time,  only  that  portion  of  the  code  was  paralle¬ 
lized.  The  inflow/outflow  particle  management,  computation  of  means,  and  other  overhead  functions  were 
performed  redundantly  on  each  processor.  The  storage  requirements  of  the  code  are  modest  enough  to 
allow  the  global  data  to  be  replicated  on  each  processor.  The  parallel  implementation  is  targeted  at  distrib¬ 
uted  memory  MIMD  (multiple  instruction,  multiple  data  set)  computers  such  as  the  Intel  iPSC/860  and  Intel 
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Delta.  The  code  was  written  in  FORTRAN  77  with  Intel’s  message  passing  calls  to  support  the  required  inter¬ 
processor  communication. 

Particles  can  be  allocated  to  processors  either  statically  or  dynamically.  In  a  static  allocation,  the 
number  of  particles  in  the  simulation  is  divided  by  the  number  of  available  processors  as  evenly  as  possible. 
In  a  dynamic  allocation,  the  particles  are  dealt  dynamically  to  processors  by  some  control  processor.  As  a 
processor  finishes  the  calculations  for  a  particular  particle,  it  is  assigned  the  job  of  solving  the  equations  for 
the  next  particle.  Note  that  dynamic  allocation,  in  which  different  processors  perform  different  tasks,  re¬ 
quires  a  MIMD  architecture.  SIMD  computers  (SIMD  stands  for  single  instruction,  multiple  data  set)  in  which 
each  processor  does  essentially  the  same  thing  at  the  same  time,  are  only  suitable  for  static  allocation. 

An  important  characteristic  of  the  algorithm  considered  here  is  that  the  solution  time  varies  signifi¬ 
cantly  from  particle  to  particle.  The  ordinary  differential  equations  for  each  particle  are  solved  by  an  explicit 
time-marching  Gear’s  method  which  adaptively  adjusts  the  step  size  to  maintain  stability.  The  solution 
times  vary  because  of  the  varying  stiffness  of  the  equations  for  different  particles.  As  a  result,  static  alloca¬ 
tion  of  equal  numbers  of  particles  to  each  processor  leads  to  poor  load  balance.  To  illustrate  this  in  specific 
terms,  we  repeat  combustion  of  CO/H2  fuel  as  represented  by  a  scheme  of  18  species  and  43  reactions 
[Correa  (1993)]  with  parameters  set  such  that  all  runs  involved  360  global  time  steps  and  762  partioles. 
Figure  10  shows  typical  results  of  static  allocation  with  10,  20,  30,  and  40  particles  assigned  per  processor. 
Notice  that  the  load  balance  does  not  improve  substantially  as  the  number  of  particles  per  processor  is  in¬ 
creased,  as  might  be  expected  from  the  law  of  averages.  In  fact,  the  imbalance  with  30  particles  per  proces¬ 
sor  is  seen  to  be  significantly  worse  than  with  either  10  or  20  particles  per  processor. 

Dynamic  allocation  through  the  use  of  a  “manager-worker”  scheme  is  found  to  lead  to  significantly 
better  load  balance,  and  to  higher  parallel  efficiencies  (i.e.,  speedup  relative  to  a  single  processor  normal¬ 
ized  by  the  number  of  processors).  One  processor  is  designated  the  “manager,”  and  given  the  job  of  deal¬ 
ing  out  the  particles  to  the  other  processors  (“workers”).  As  a  worker  processor  finishes  calculating  the  time 
step  for  its  particle,  it  signals  the  manager  that  it  is  available  to  perform  another  calculation.  In  this  manner,  a 
processor  that  is  assigned  particles  that  can  be  solved  quickly  will  calculate  more  particles,  while  a  proces¬ 
sor  that  is  assigned  particles  that  take  a  long  time  to  solve  will  calculate  fewer  particles.  Figure  1 1  shows  the 
number  of  particles  dynamically  assigned  to  each  processor  for  the  same  test  case  considered  in  Fig.  10. 
The  number  of  particles  assigned  is  not  at  all  constant,  but  is  seen  to  vary  significantly  from  processor  to 
processor.  When  the  average  number  of  particles  per  processor  exceeds  ten  or  so,  dynamic  allocation 
smooths  out  the  load  balance  well  (Fig.  1 2) .  The  maximum  time  required  by  any  processor  is  also  seen  to  be 
significantly  smaller  with  dynamic  allocation  than  the  maximum  time  required  by  any  processor  with  static 
allocation,  as  comparison  of  Figs.  10  and  12  shows.  The  comparison  is  striking  in  the  case  with  480  par¬ 
ticles,  where  dynamic  allocation  reduces  the  maximum  time  from  22.5  seconds  to  about  19  seconds. 

For  dynamic  allocation  to  be  efficient,  each  processor  must  have  equal  access  to  the  data  needed  to 
perform  the  requested  calculation.  The  algorithm  under  consideration  here  requires  relatively  little  total  stor¬ 
age,  allowing  the  global  data  to  be  replicated  on  each  processor  and  hence  making  dynamic  allocation 
feasible.  One  drawback  of  the  Intel  iPSC/860  is  that  only  a  single  process  can  run  on  each  processor.  As  a 
result,  the  manager  processor  sits  idle  much  of  the  time  while  the  worker  processors  perform  their  computa¬ 
tions,  because  it  is  not  possible  to  start  an  additional  worker  process  on  the  manager  processor.  A  practical 
consequence  of  this  is  that  on  a  computer  with  Np  processors,  only  Np-1  processors  are  solving  the  par- 
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tide  equations.  This  is  clearly  more  of  a  drawback  on  small  computers  than  on  large  massively  parallel  ma¬ 
chines.  The  improved  load  balancing  achieved  by  dynamic  allocation  more  than  offsets  this  disadvantage 
on  the  larger  machines.  Dynamic  allocation  is  superior  to  static  allocation  when  six  or  more  processors  are 
used.  The  falloff  in  the  improvement  achieved  on  the  Delta  as  the  number  of  processors  is  increased  is  a 
result  of  the  number  of  particles  per  processor  becoming  small  (here,  <  10),  which  reduces  the  effective¬ 
ness  of  dynamic  allocation  in  smoothing  out  the  load  imbalance  between  processors  as  discussed  above. 

Another  important  consideration  in  the  parallelization  of  Monte-Carlo  algorithms  is  the  generation  of 
random  numbers  on  a  parallel  computer  versus  a  serial  computer.  One  of  the  objectives  in  parallelizing  the 
PaSR  model  was  to  ensure  that  the  both  the  convergence  path  and  the  final  solution  were  independent  of 
the  number  of  processors  used.  This  requires  that  the  same  particles  are  removed  from  the  reactor  as  new 
particles  are  introduced  on  each  and  every  processor  at  every  stage  of  the  calculation.  This  was  accom¬ 
plished  in  an  almost  trivial  manner  here  by  redundantly  performing  the  particle  removal  on  each  processor, 
and  by  starting  the  same  random  number  generator  on  each  processor  with  the  same  seed  so  that  identical 
sequences  of  random  numbers  were  generated.  Had  this  portion  of  the  calculation  also  been  performed  in 
parallel,  special  care  would  have  been  required  to  guarantee  identical  results. 

The  calculation  of  the  mean  at  the  end  of  the  time  step  requires  global  communication  between  the 
processors  since  each  processor  only  updates  the  solution  at  the  time  step  of  interest  for  a  sub-ensemble 
of  particles.  There  are  two  different  ways  of  calculating  the  mean  that  can  be  used.  One  method  is  to  have 
each  processor  compute  the  sum  of  all  of  the  properties  on  its  sub-ensemble,  and  then  globally  exchange 
the  sums.  Another  method  is  to  have  each  processor  perform  a  global  shuffle  of  each  of  the  columns  of  the 
solution  vector  it  computed  with  each  of  the  other  processors,  so  that  each  processor  ultimately  recon¬ 
structs  the  entire  solution  vector.  At  this  point,  the  mean  properties  can  all  be  computed  redundantly  by 
each  processor,  using  the  same  code  as  in  the  serial  implementation.  The  first  method  is  more  computation¬ 
ally  efficient  because  the  computation  of  the  means  is  done  partially  in  parallel,  and  because  less  data  must 
be  communicated.  However,  the  method  suffers  the  potential  drawback  of  roundoff  error  because  the  sums 
on  the  different  processors  are  not  all  computed  in  the  same  order  and  the  results  may  be  slightly  different. 
The  second  method  requires  more  comunication,  but  ensures  identical  results  on  each  processor. 

Timing  results  indicate  that  the  parallel  inefficiency  observed  is  caused  primarily  by  load  imbalance, 
and  that  the  time  required  to  globally  communicate  the  results  and  compute  the  means  is  insignificant.  Con¬ 
sequently,  either  method  of  computing  the  means  appears  adequate  for  this  application.  The  second  meth¬ 
od  was  used  here  because  it  required  fewer  modifications  to  the  code  and  ensured  identical  results. 

To  verify  the  parallelization  and  demonstrate  the  speedup,  a  CO/H2  scheme  per  Table  I  was  calcu¬ 
lated  on  (i)  four  serial  computers:  a  Sun  Sparcstation  1 ,  a  Convex  C21 0,  a  Cray  Y-MP,  and  a  Cray  Y-MP 
C90  (hereafter  called  a  Cray  C90),  and  (ii)  two  parallel  computers:  a  hypercube-based  16-node  Intel 
iPSC/860,  and  the  mesh-based  512-node  Intel  Touchstone  Delta  prototype.  The  Cray  runs  were  made 
with  the  original  version  of  the  code,  and  repeated  with  aversion  optimized  in  a  machine -specific  manner. 

Examining  the  runtimes  (cpu  time),  the  PaSR  code  on  the  16-node  Intel  hypercube  was  found  to  run 
about  1.4  times  faster  than  on  the  Cray  C90  and  about  17.1  times  faster  than  on  the  CONVEX.  The  run  on 
256  processors  of  the  Intel  Delta  took  26  minutes.  This  result  was  17.7  times  faster  than  the  run  of  the  original 
code  on  the  C90,  and  10.7  times  faster  than  the  run  of  the  modified  code  on  the  C90.  Although  a  single 
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processor  of  the  Intel  computers  is  about  a  factor  of  ten  slower  than  the  Cray  C90,  the  PaSR  algorithm  does 
not  “vectorize”  well  and  therefore  the  principal  feature  of  the  latter  machine  cannot  be  exploited.  It  should 
be  mentioned  that  this  is  not  characteristic  of  all  classes  of  equations;  for  example,  conventional  computa¬ 
tional  fluid  dynamics  (CFD)  codes  are  significantly  accelerated  on  vector  computers. 

Summarizing  the  results  of  parallelization.  Fig.  13  indicates  that  the  rate  of  computation  (“Mflop”  = 
millions  of  floating  point  operations  per  second)  is  linear  in  the  number  of  processors,  up  to  the  level  of  64 
processors  on  the  Intel  machines.  A  yield  of  approximately  4  Mflop/processor  is  obtained  on  the  iPSC/860 
processors  in  this  range.  The  speedup  is  less  than  linear  as  the  number  of  particles  per  processor  becomes 
less  than  10,  because  then  the  dynamic  allocation  is  less  able  to  achieve  good  load-balancing.  The  initial 
conditions  on  the  particle  in  each  global  time  step  greatly  affect  the  stiffness  of  the  governing  o.d.e.’s,  and 
so  the  computational  time.  With  a  large  number  of  particles  per  processor,  the  total  computational  load  be¬ 
comes  similar  for  each  processor;  at  the  other  extreme,  the  load  can  be  very  disparate.  In  this  limit,  most  of 
the  processors  may  have  to  “wait”  for  the  few  that  are  integrating  the  stiffest  sub-sets  of  the  equations.  On 
the  other  hand,  had  about  5,000  particles  been  used  in  Fig.  13  (note  that  115,000  particles  were  used  in  the 
combined  CFD/pdf  transport  calculation  of  Section  2.2,  albeit  with  simple  chemistry)  a  rate  of  about  2 
GFIops  (4  Mflops  x  512  processors  of  the  Delta  machine)  could  have  been  attained. 

It  is  clear  that  distributed  memory  MIMD  parallel  architectures  are  well  suited  to  the  particle -tracking 
PaSR  algorithm.  This  remains  true  so  long  as  the  turbulence  model  does  not  require  continuous  commu¬ 
nication  between  “particles.”  Thus,  the  lEM  model  is  preferable  to  pair-exchange  models  on  the  parallel 
computers. 

2.6  CH4/Air  System 

We  have  used  the  PaSR  model  to  study  a  27-species/77- reaction  methane  oxidation  scheme  per 
Table  II  [Correa  and  Braaten  (1993)].  The  calculation  involved  600  to  700  particles  and  therefore  —20,000 
o.d.e.’s.  Figure  14  shows  typical  results  for  a  30  atm,  1200K  (inlet)  premixed  methane-air  system.  Figure 
14  (a)  shows  the  convergence  of  ensemble-mean  Yqo  to  a  stochastic  steady-state;  the  degree  of  turbu¬ 
lence  in  the  time -histories  of  species  is  greatest  for  CO  (among  major  species)  as  blowout  is  approached. 
Figure  1 4(b)  shows  that  the  PaSR  solutions  properly  meet  the  limits  imposed  by  PSR  and  PFR  models;  note 
that  with  this  choice  of  inlet  conditions,  reactor  volume  and  mass  flow  (such  that  the  residence  time  t  turns 
out  to  be  approximately  2  ms),  the  flame  blows  out  in  the  low  frequency  limit.  Figure  14  (c)  shows  that  the 
“skeletal”  mechanism  (used  by  Seshadri  and  Peters  (1990)  as  the  basis  for  the  4-step  reduced  scheme 
in  laminar  flames)  cannot  reproduce  the  full  scheme's  results  at  low  frequency,  where  ignition  chemistry 
becomes  important.  Figure  14  (d)  shows  atypical  scatterplot .  here  for  NOx  although  anyone  of  the  27  spe¬ 
cies  or  temperature  could  have  been  shown.  The  scatterplots  can  be  used  to  derive  pdf’s,  correlations, 
rms’s,  and  other  quantities. 

We  have  only  briefly  summarized  the  development  and  some  key  conclusions  of  the  PaSR  model. 
Many  more  results  and  sensitivity  analyses  are  given  in  Correa  and  Braaten  (1993),  Correa  (1993),  and  Cor¬ 
rea  (1994  a).  The  PaSR  model  has  also  been  adopted  as  an  investigative  tool  by  other  researchers,  e.g, 
at  the  University  of  Sydney  and  at  UC-Berkeley. 


Table  I.  Kinetic  Mechanism  Used  for  C0/H2-Air  Combustion. 
(Forward  rate  constant  kf  =  A  moles-cm-s-K;  E  units:  cal/mole) 


No.  Reaction 

A 

b 

E 

1.  C0+0+M=C02+M 

3.20E+13 

0.0 

-4200.0 

2.  CO+OH=C02+H 

1.51E+07 

1.3 

-758.0 

3.  C0+02=C02+0 

1.60E+13 

0.0 

41000.0 

4.  H02+C0=C02  +  0H 

5.80E+13 

0.0 

22934.0 

5.  H2+02=20H 

1.70E+13 

0.0 

47780.0 

6.  0H+H2=H20+H 

1.17E+09 

1.3 

3626.0 

7.  H+02=0H+0 

2.00E+14 

0.0 

16800.0 

8.  0+H2=OH  +  H 

1.80E  +  10 

1.0 

8826.0 

9.  H+02  +  M  =  H02  +  M 

2.10E+18 

-1.0 

0,0 

H20  enhanced  by  21 .0,  C02  enhanced  by  5.0,  H2  enhanced  by  3.3, 
CO  enhanced  by  2.0,  02  enhanced  by  0.0,  N2  enhanced  by  0.0 

10.  H+02+02=H02+02 

6.70E+19 

-1.4 

0.0 

11.  H+02+N2  =  H02  +  N2 

6.70E+19 

-1.4 

0.0 

12.  0H  +  H02  =  H20  +  02 

5.00E+13 

0.0 

1000.0 

13.  H  +  H02=20H 

2.50E  +  14 

0.0 

1900.0 

14.  0  +  H02=02  +  0H 

4.80E  +  13 

0.0 

1000.0 

15.  20H  =  0  +  H20 

6.00E  +  08 

1.3 

0.0 

16.  H2+M  =  H+H  +  M 

2.23E  +  12 

0.5 

92600.0 

H20  Enhanced  by  6.0,  H  enhanced  by  2.0,  H2  enhanced  by  3.0 

17.  02+M=0  +  0  +  M 

1.85E  +  11 

0.5 

95560.0 

18.  H+0H+M  =  H20+M 

7.50E+23 

-2.6 

0.0 

H20  enhanced  by  21.0 

19.  H+H02=H2+02 

2.50E+13 

0.0 

700.0 

20.  H02  +  H02  =  H202+02 

2.00E  +  12 

0.0 

0.0 

21.  H202+M=0H+0H  +  M 

1.30E+17 

0.0 

45500.0 

22,  H202+H  =  H02+H2 

1.60E+12 

0.0 

3800.0 

1.00E4-13 


0.0 


1800.0 


CPU  TIME  (sec  /  time  step) 


PROCESSOR  NUMBER 


Figure  10.  Load  Imbalance  resulting  from  static  allocation  of  particles  to  processors. 
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PROCESSOR  NUMBER 


Figure  12,  Load  balance  resulting  from  dynamic  allocation  per  Fig.  11. 
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M  FLOPS 


4  MFLOPS  /  PROCESSOR 
REFERENCE  LINE 


SERIAL  MACHINES 
CRAY  C90 
CRAYYMP 
C210 
SPARC  1 


PARALLEL  MACHINES 

■  iPSC  860  (16  proc.) 
H  DELTA  (64  proc.) 
A  DELTA  (128  proc.) 
X  DELTA  (256  proc.) 


NUMBER  OF  PROCESSORS 


Figure  13,  Performance  of  PaSR  algorithm  on  serial  and  parallel  computers  for  CO/H2  fuel,  1  atm,  1000K  inlet,  762  particles 


Table  II.  “Full”  Kinetic  Scheme  used  for  Lean  Methane  Combustion  in  PaSR  Model. 


(Forward  rate  constant  kf  =  A  moles-cm-s-K;  E  units:  cal/mole) 


No. 

Reaction 

A 

b 

E 

1. 

CH4(+M)=CH3+H(+M) 

6.30E+14 

0.0 

104000.0 

Low  pressure  limit: 

6.30E-03 

0.0 

18000.0 

2. 

CH4+02=CH3+H02 

7.90E+13 

0.0 

56000.0 

3. 

CH4+H=CH3+H2 

2.20E+04 

3.0 

8750.0 

4. 

CH4+0=CH3+0H 

1.60E+06 

2.4 

7400.0 

5. 

CH4+0H=CH3+H20 

1.60E+06 

2.1 

2460.0 

6. 

CH20+0H=HC0+H20 

7.53E+12 

0.0 

167.0 

7. 

CH20+H=HC0+H2 

3.31E+14 

0.0 

10500.0 

8. 

CH20+M=HC0+H+M 

3.31E+16 

0.0 

81000.0 

9. 

CH20+0=HC0+0H 

1.81E+13 

0.0 

3082.0 

10. 

HC0+0H=C0+H20 

5.00E+12 

0.0 

0.0 

11. 

HCO+M=H+CO+M 

1.60E+14 

0.0 

14700.0 

12. 

HCO+H==CO+H2 

4.00E+13 

0.0 

0.0 

13. 

HC0+0=0H+C0 

1.00E+13 

0.0 

0.0 

14. 

HC0+02=H02+C0 

3.00E+12 

0.0 

0.0 

15. 

C0+0+M=C02+M 

3.20E+13 

0.0 

-4200.0 

16. 

C0+0H=C02+H 

1.51E+07 

1.3 

-758.0 

17. 

C0+02=C02+0 

1.60E+13 

0.0 

41000.0 

18. 

CH3+02-CH30+0 

7.00E+12 

0.0 

25652.0 

19. 

CH30+M=CH20+H+M 

2.40E+13 

0.0 

28812.0 

20. 

CH30+H=CH20+H2 

2.00E+13 

0.0 

0.0 

21. 

CH30+0H=CH20+H20 

1.00E+13 

0.0 

0.0 

22. 

CH30+0=CH20+0H 

1.00E+13 

0.0 

0.0 

23. 

CH30+02=CH20+H02 

6.30E+10 

0.0 

2600.0 

24. 

CH3+02=CH20+0H 

5.20E+13 

0.0 

34574.0 

25. 

CH3+0=CH20+H 

6.80E+13 

0.0 

0.0 

26. 

CH3+0H  =  CH20+H2 

7.50E+12 

0.0 

0.0 

27. 

CH2+H=CH+H2 

4.00E+13 

0.0 

0.0 

28. 

CH2+0=C0+H+H 

5.00E+13 

0.0 

0.0 

29. 

CH2+02=C02+H+H 

1.30E+13 

0.0 

1500.0 

30. 

CH+0=C0+H 

4.00E+13 

0.0 

0.0 

31. 

CH+02=C0+0H 

2.00E+13 

0.0 

0.0 
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32. 

CH2CO  +  H=CH3+CO 

7.00E+12 

0.0 

3000.0 

33. 

CH2C0+0-HC0+HC0 

2.00E+13 

0.0 

2300.0 

34. 

CH2C0+0H=CH20+HC0 

1.00E+13 

0.0 

0.0 

35. 

CH2CO  +  M=CH2+CO+M 

1.00E+16 

0.0 

59250.0 

36. 

HCCO  +  H=CH2+CO 

3.00E+13 

0.0 

0.0 

37. 

HCC0+0=C0+C0+H 

1.20E+12 

0.0 

0.0 

38. 

H02+C0=C02+0H 

5.80E+13 

0.0 

22934.0 

39. 

H2+02=20H 

1.70E+13 

0.0 

47780.0 

40. 

0H+H2=H20+H 

1.17E+09 

1.3 

3626.0 

41. 

H+02=0H+0 

2.00E  +  14 

0.0 

16800.0 

42. 

0+H2=0H+H 

1.80E+10 

1.0 

8826.0 

43. 

H+02+M=H02+M 

2.10E  +  18 

-1.0 

0.0 

H20  Enhanced  by  2.1 ,  C02 

Enhanced  by  5.0,  H2  Enhanced  by  3.3, 

CO  Enhanced  by  2.0,  02  Enhanced  by  0.0,  N2  Enhanced  by  0.0 

44. 

H+02+02-H02+02 

6.70E+19 

-1.4 

0.0 

45. 

H+02+N2=H02+N2 

6.70E  +  19 

-1.4 

0.0 

46. 

0H  +  H02=H20+02 

5.00E+13 

0.0 

1000.0 

47. 

H+H02=20H 

2.50E+14 

0.0 

1900.0 

48. 

0+H02=02+0H 

4.80E+13 

0.0 

1000.0 

49. 

20H=0+H20 

6.00E+08 

1.3 

0.0 

50. 

H2+M=H+H+M 

2.23E+12 

0.5 

92600.0 

H20  Enhanced  by  6.0,  H  Enhanced  by  2.0,  H2  Enhanced  by  3.0 

51. 

02+M=0+0+M 

1.85E+11 

0.5 

95560.0 

52. 

H+0H+M-H20+M 

7.50E+23 

-2.6 

0.0 

H20  Enhanced  by  2.0 

53. 

H+H02=H2+02 

2.50E+13 

0.0 

700.0 

54. 

H02+H02=H202+02 

2.00E  +  12 

0.0 

0.0 

55. 

H202+M=0H+0H+M 

1.30E+17 

0.0 

45500.0 

56. 

H202+H=H02+H2 

1.60E  +  12 

0.0 

3800.0 

57. 

H202+0H=H20+H02 

1.00E+13 

0.0 

1800.0 

58. 

N+02=N0+0 

6.40E+09 

1.0 

6280.0 

59. 

N+OH=NO+H 

3.80E  +  13 

0.0 

0.0 

60. 

N+N0=N2+0 

3.30E+12 

0.3 

0.0 

61. 

N+C02=N0+C0 

1.90E  +  11 

0.0 

3400.0 

62. 

N0+H02=N02+0H 

2.10E+12 

0.0 

-480.0 

63. 

N02+M=N0+0+M 

1.10E+16 

0.0 

66000.0 

64. 

N02+H=:N0+0H 

3.50E+14 

0.0 

1500.0 

65. 

N02+0=N0+02 

1.00E+13 

0.0 

600.0 
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66.  HNO  +  M  =  H  +  NO+M  1.50E+16  0.0  48680.0 

H20  Enhanced  by  6.0,  H2  Enhanced  by  2.0,  02  Enhanced  by  2.0, 


67. 

N2  Enhanced  by  2.0 

HNO+H=H2+NO 

68. 

HN0+0H=N0+H20 

69. 

N20=N2+0 

70. 

N20+H=N2+0H 

71. 

N20+0=N0+N0 

72. 

N20+0=N2+02 

73. 

NCO+M=N+CO+M 

74. 

NC0+0=N0+C0 

75. 

NCO+OH=NO+CO+H 

76. 

NCO+N=N2+CO 

77. 

NC0+N0=N20+C0 

5.00E+12 

0.0 

0.0 

3.60E+13 

0.0 

0.0 

2.82E+16 

-1.6 

62130.0 

7.60E+13 

0.0 

15200.0 

1.00E+14 

0.0 

28200.0 

1.00E+14 

0.0 

28200.0 

3.10E+16 

-0.5 

48000.0 

5.60E+13 

0,0 

0.0 

1.00E+13 

0.0 

0.0 

2.00E+13 

0.0 

0.0 

1.00E+13 

0.0 

-390.0 

MEAN  MASS  FRACTION  OF  CO,  ppm 


(a)  Convergence  of  PaSR  to  stochastic  steady  state. 

Notice  the  greater  level  of  turbulence  in  CO  (and  others) 
at  the  lower  frequency  (nearer  blowout  at  these  conditions). 


(b)  Mean  PaSR  temperature  and  limit  values.  Gas-turbine 
combustors  are  being  pushed  towards  the  higher  mixing  rates, 
or  else  the  blowout  (extinction)  threshold  will  be  crossed. 
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(c)  Comparison  of  skeletal  and  full  schemes.  The  25-step  (d)  Scatterplot  of  NO  in  the  1000  Hz  ensemble  after  stochastic 
skeletal  scheme  blows  out  prematurely.  The  derivative  steady  state  is  reached.  Note  wide  distribution  of  NO.  Correspond- 
4 -step  reduced  scheme  blew  out  at  all  frequencies.  ing  scatterplots  show  wide  variations  in  CO,  CHx  and  temperature,  i 


Figure  14.  Selected  results  from  the  "PaSR”  model. 


28 


2.7  Assessment  Of  A  3-Variable  Reduced  Scheme  For  Complex  Fuel 


Most  research  in  modeling  is  conducted  with  simple  fuels,  such  as  CO/H2  mixtures  or  CH4.  However, 
jet  engines  burn  higher  hydrocarbons  of  typical  composition  CnHm  with  n  =  10  and  m==20.  The  detailed 
kinetics  of  the  oxidation  mechanism  are  unknown.  Even  if  they  were  known,  computationally  tractable  kinet¬ 
ic  schemes  would  be  required  for  these  fuels.  Quasi-global  schemes  have  been  used  in  the  past  for  engine 
applications  [Edelman  and  Harsha  (1978)],  as  well  as  for  stirred  reactors  [Duterque  (1981)]  and  laminar 
flames  [Westbrook  and  Dryer(1981)].  Recent  experimental  data  offer  much  detail  on  the  intermediates 
formed  in  oxidation  of  kerosene  in  a  jet-stirred  reactor,  thereby  guiding  the  formulation  of  pyrolysis  models 
[Gueret  et  al.  (1990)]. 

The  four-  and  five -variable  schemes  [e.g.,  Seshadri  and  Peters  (1990)]  are  not  worth  the  additional 
computational  burden,  given  the  empiricism  in  the  global  pyrolysis  rate  for  complex  fuels  as  well  as  the 
shortcomings  of  the  four-step  scheme  demonstrated  above  in  the  PaSR.  Hence,  here  we  will  develop  a 
three-variable  scheme  that  explicitly  permits  and  retains  a  global  pyrolysis  rate. 

To  realistically  assess  a  reduced  kinetic  scheme,  the  testbed  must  reproduce  the  microscale  turbu¬ 
lent  environment  of  the  burner,  or  at  least  that  of  the  CFD  methodolgy  in  which  the  reduced  scheme  will 
be  employed.  To  this  end,  the  Partially  Stirred  Reactor  ("PaSR”)  model  is  used,  taking  advantage  of  charac¬ 
teristics  such  as  (i)  its  behavior  is  bounded  by  the  perfectly  stirred  reactor  (PSR)  and  the  plug  flow  reactor 
(PFR),  providing  a  set  of  checks  on  the  model;  (ii)  at  the  high-frequency  end,  it  provides  a  turbulence  envi¬ 
ronment  typical  of  highly  turbulent  combustors  [Correa  (1989)]  namely,  those  in  the  distributed  reaction 
zone  regime  rather  than  in  the  v\ieaker  laminar  flamelet  regime,  and  (iii)  the  joint  velocity— scalar(s)  pdf  trans¬ 
port  equation  for  multi-dimensional  Monte  Carlo  CFD  degenerates  to  the  PaSR  within  each  computational 
cell. 


Following  assessment  in  the  PaSR  over  the  entire  range  of  mixing  frequencies  possible,  a  reduced 
scheme  can  be  implemented  in  the  multi-dimensional  pdf/CFD  model,  and  those  predictions  can  be 
compared  with  data  (from  “real— world”  inhomogeneous  flow).  This  approach  will  provide  the  combustor 
designer  with  a  predictive  methodology  that  has  been  validated,  to  the  maximum  extent,  at  each  step. 

Since  detailed  kinetic  schemes  for  complex  fuels  are  unavailable,  as  discussed  above,  the  baseline 
or  “starting”  kinetic  scheme  is  initiated  by  an  irreversible  global  pyrolysis  reaction  [Correa  (1994  b)] 

CnHm  +  (n/2)  O2  n  CO  +  (m/2)  H2 
whose  rate  can  be  a  function  of  the  equivalence  ratio;  chain -branching 
CO  +  OH  ^  CO2  +  H 
H2  +  0  OH  +  H 
H  +  O2  ^  OH  +  0 
OH  +  OH  H2O  +  0 
and  recombination  reactions  such  as 
H2  +  M->H  +  H  +  M 
O2  +  M— »  0  +  0  +  M 
H  +  OH  +  M  ^  H2O  +  M 

Standard  rates  are  used  for  the  elementary  reactions  (5)-(11). 


(4) 

(or  “shuffle”)  reactions  such  as 

(5) 

(6) 

(7) 

(8) 

(9) 

(10) 
(11) 
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Multi-step  schemes  contain  too  many  variables  (species  plus  enthalpy)  to  be  tractable  in  CFD. 
Hence  the  number  of  variables  must  be  reduced.  Since  the  multi-step  pyrolysis  of  lower  hydrocarbons 
has,  of  necessity,  already  been  replaced  by  the  assumed  single-step  pyrolysis  of  the  complex  fuel,  a  good 
start  has  been  made  toward  a  reduced  mechanism.  Here  partial  equilibrium  in  the  radical  pool  is  used  to 
further  eliminate  (chain-branching)  steps. 

For  purposes  of  this  development,  we  assume  that  “air  ”  consists  of  21  %  O2  and  79%  N2,  by  volume, 
and  that  the  initial  reactants  are  exclusively  the  fuel  (CnHm)  and  air.  Furthermore,  the  only  species  in  the 
system  are  CnH^,  O2,  N2,  CO,  H2,  CO2,  H2O,  0,  OH,  and  H.  Define  Wf=(nWc+mWH)  and 
Wa=  (W02+3.76W1S12),  where  Wj  is  the  molecular  weight  of  species  “i”  (except  that  Wg  is  not  the  mean  mo¬ 
lecular  weight  of  air).  This  system  is  described  in  terms  of  three  variables:  the  mixture  fraction  the  fuel 
mass  fraction  Yf,  and  a  composite  radical  pool  variable  Y*.  The  ranges  of  the  three  variables  Yf  and  Y* 
follow  and  will  show  that  all  possible  thermochemical  states  fall  within  a  tetrahedron  in  this  three-dimen¬ 
sional  (|-Yf-Y*)  space. 

The  first  variable  is  a  conserved  scalar  called  the  mixture  fraction  which  is  derived  by  normalizing 
the  elemental  mass  fractions,  by  the  values  in  the  fuel  and  air  streams. 

ZrZ’ 


Superscripts  a  and  f  refer  to  the  air  and  fuel  streams,  respectively.  ^  is  conserved  because  elements  are 
unchanged  by  chemical  reaction;  likewise,  the  total  (chemical  plus  sensible)  enthalpy  is  conserved  under 
reaction  and  can  be  mapped  into 


By  construction,  therefore.  0  <  ^  <  1, 

The  second  variable  is  the  fuel  mass  fraction  Yf  <  Yf  <  Yf^^^(^).  The  upper  bound  is  set  by 

having  all  the  elements  present  in  the  initial  reactants,  i.e,,  Yf^^^(§)  =  The  lower  bound  Yf^‘'^(5)  is  set  by 
the  availability  of  O2  to  convert  the  maximum  possible  amount  of  C  to  CO  and  H2  (not  to  CO2  and  H2O). 

The  third  variable  is  a  composite  Y*,  defined  such  that  o*  =  oh2  +  where  Oj  =  Yj  /  Wj  is  the  mole 
number  of  species  “j.”  The  bounds  on  Y*  follow  from  the  bounds  on  CO  and  H2  [Correa  (1993  c)]. 

By  considering  the  bounds  on  these  three  variables,  the  allowable  ^-Yf-Y*  space  within  which  the 
reduced  chemistry  scheme  must  lie  is  obtained  as  a  tetrahedron  (Fig.  15). 

Post- steady  “State  ensemble-mean  quantities  obtained  from  the  PaSR  model,  using  the  starting  and 
the  reduced  schemes,  are  shown  in  Figure  1 6.  In  each  case,  the  independently  computed  high-frequency 
limit  (PSR)  and  low-frequency  limit  (PFR  convoluted  with  the  theoretical  age  distribution)  are  also  shown. 
All  PaSR  predictions  lie  smoothly  between  these  limits,  providing  a  set  of  checks  on  the  model.  Several 
points  can  be  made: 

(1)  The  global  pyrolysis  rate  In  the  starting  scheme  does  not  yield  complete  combustion,  appropriate  for 
a  heavy  hydrocarbon  under  present  circumstances.  The  low-frequency  limit  has  only  a  360K  temperature 
rise;  mixing  is  required  to  accelerate  ignition  of  the  incoming  reactants,  but  even  in  the  high-frequency 
mixing  limit,  the  PSR  temperature  rise  is  only  860K  and  not  the  1500K  equilibrium  rise.  The  assumed  rate 
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can  be  altered  to  produce  other  results,  if  desired. 

(2)  The  agreement  between  the  starting  scheme  and  the  reduced  scheme  is  excellent  on  the  mean  temper¬ 
ature  (Fig.  16a)  and  on  the  mean  fuel  mass  fraction  (Fig.  16b),  at  ali  frequencies  from  10  to  1000  Hz.  The 
agreement  on  fuel  mass  fraction  is  not  surprising  since  fuel  pyrolysis  is  explicitiy  recognized  as  a  degree- 
of-freedom  in  the  reduced  scheme.  The  latter  is  a  desirable  feature. 

(3)  The  agreement  on  mean  Y*,  which  is  the  combined  variable  based  on  the  CO  and  H2  mole  numbers, 
is  not  as  good  (Fig.  16c).  Potential  contributors  to  this  discrepancy  are  the  integration  error,  the  coarseness 
in  the  look-up  table,  and  the  assumption  of  partial  equilibrium.  These  factors  can  be  examined  in  turn.  First, 
the  integration  time  step  ftt  is  small  enough  (8t  =0.2  ps  « 1  /  w^^ax)  to  resolve  the  fastest  rates  in  the  system; 
in  fact,  simulations  with  smaller  time  steps  gave  the  same  results  to  several  significant  figures.  Second,  ex¬ 
amination  of  the  1000  Hz  ensemble  shows  that  most  of  the  particles  are  adjacent  to  the  minimum  Y*  bound¬ 
ary;  hence,  the  coarseness  of  the  grid  in  the  Y*  direction  (21  uniformly  spaced  points  between  Y*-'^''^  and 
Y*,max)  contributes  to  the  discrepancy.  This  error  can  be  reduced  at  the  straightforward  expense  of  adding 
points  to  the  table  in  the  Y*  direction.  Third,  the  assumption  of  partial  equilibrium  in  reactions  5-8  is  not 
strictly  correct  under  conditions  where  significant  amounts  of  fuel  are  present;  hence,  it  is  responsible  in 
part  for  the  discrepancies  between  the  starting  and  reduced  schemes. 

On  balance,  the  reduced  scheme  exhibits  encouraging  performance.  Recent  data  on  the  intermediates 
formed  in  the  oxidation  of  kerosene  in  a  jet— stirred  reactor  [Gueret  et  al.  (1990)]  will  guide  further  develop¬ 
ment. 


Figure  15.  Depiction  of  ailowable  |-Yf-o*  space  in  |-Yf  plane. 
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ENSEMBLE  MEAN  TEMPERATURE,  K 


z 


TURBULENCE  FREQUENCY  Hz 


(a)  Mean  temperature. 


(b)  Mean  fuel  mass  fraction. 


(c)  Mean  Y*,  combined  variable  based  on  the  CO  and  H2  mole  numbers. 

Figure  16.  Comparison  of  ensemble-mean  quantities  obtained  from  the  PaSR  model  using  the  starting  and 
the  reduced  schemes  for  a  C^Hm  fuel.  The  PFR  and  PSR  limits  were  obtained  using  the  starting  scheme 
and  independent  codes. 
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2.8  Direct  Comparison  of  Mixing  Models 

Because  reactions  introduce  at  least  one  chemical  time  scale  into  the  flow,  the  ratio  of  a  characteristic 
mixing  time  to  a  characteristic  chemical  time  (the  Damkohler  number,  Da)  becomes  useful  to  classify  the 
regime.  If  Da>1  the  system  is  in  the  flamelet  regime  where  mixing  models  can  be  very  inaccurate  [Chen 
and  Kollmann  (1990)].  Concerns  arise  because  pair-exchange  models  treat  mixing  and  chemistry  as  se¬ 
quential  rather  than  simultaneous  events,  which  can  result  in  non-physical  behavior.  For  example,  consid¬ 
er  two  fluid  particles  which  are  outside  the  scalar(s)’  range  of  chemical  reactivity,  both  individually  and  in 
the  mean.  Some  chemical  reactions  should  occur  in  the  process  of  mixing  through  the  reactive  regime; 
however,  some  models  will  return  the  non -reactive  mean  and  hence  mixing  will  take  place  without  chemi¬ 
cal  reaction.  If  Dacl  the  system  is  in  the  “distributed  combustion”  regime,  and  reactions  do  not  influence 
the  turbulence  to  a  significant  extent.  In  reality,  combustion  proceeds  through  a  multi-step  reaction  mecha¬ 
nism  with  a  wide  range  of  Damkohler  numbers  potentially  spanning  the  gamut  from  distributed  to  flamelet 
behavior.  Hence,  little  of  useful  generality  can  be  said  regarding  scalars  which  react  according  to  a  realistic 
chemical  kinetic  scheme.  Numerical  experiments  are  needed  to  assess  mixing  models  in  this  context. 

The  lEM  model  was  advanced  as  being  appropriate  in  Lagrangian  calculations  (such  as  particle 
tracking  Monte  Carlo  calculations)  of  turbulent  reacting  flow  [Meyers  and  O’Brien  (1980)].  However,  it  has 
certain  peculiarities  such  as  the  shape-preserving  nature  of  the  relaxation  of  the  initial  pdf  of  a  conserved 
scalar  undergoing  mixing.  Hence,  it  is  of  interest  to  compare  it  with  other  mixing  models  in  the  context  of 
turbulent  combustion.  In  particular,  self- igniting  combustion  is  a  good  test  both  for  practical  reasons  — 
many  important  combustion  systems  such  as  those  in  gas-turbines  are  aerodynamically  stabilized  by  in¬ 
tense  mixing  of  products  with  reactants  —  as  well  as  theoretical  reasons,  since  a  different  set  of  time  scales 
may  be  operative  in  the  ignition  phase. 

The  objective  of  this  part  of  the  study  was  to  compare  the  Curl  (1963)  model,  a  modified  Curl  model 
developed  by  Dopazo  (1979)  and  Janicka  et  al.  (1979),  and  the  “Interaction-by-Exchange-with-the- 
Mean”  (lEM)  mixing  model  [Borghi  (1988)]  which  was  introduced  above. 

(i)  The  Curl  (1963)  model  randomly  chooses  Nmix  pairs  of  particles,  mixes  them  by  averaging  their 
scalar  values,  and  finally  returns  them  to  the  ensemble.  The  mean  is  unchanged,  whereas  the  variance  de¬ 
creases,  in  the  inert  case.  Illustrating  the  case  for  a  conserved  scalar  consider  the  mixing  of  one  pair  (out 
of  Nmix  pail's  of  particles  which  are  mixed  at  each  time  step,  where  Nmix 's  given  below).  Let  the  two  partici¬ 
pating  particles  have  the  values  51  and  52  at  beginning  of  the  time  step.  Then,  according  to  the  Curl 
model,  the  values  at  the  end  of  the  step  (superscript  “new”)  are 

tnew  ^  tnew  ^  ^  (14a) 

To  produce  the  correct  decay  rate  for  the  variance  of  a  conserved  scalar/ 

N„,  =  2  Np  0)  At  (14b) 

Weaknesses  of  the  Curl  model  have  been  described  in  detail,  in  the  case  of  the  evolution  of  the  pdf 
of  a  conserved  scalar  [Pope  (1982)].  In  particular,  the  Curl  model  never  yields  a  continuous  distribution 
when  starting  with  a  pdf  which  is  composed  of  delta  functions;  hence,  the  higher  moments  of  the  pdf  are 
never  correct. 
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Following  mixing,  the  mass  fractions  and  temperature  of  each  particle  in  the  ensemble  are  advanced 
in  time  by  integration  of  the  chemical  kinetics  equations. 


(ii)  The  Modified  Curl  model  also  randomly  chooses  pairs  of  particles,  but  mixes  them  by  averaging 
their  scalar  values  in  a  weighted  manner  [Dopazo  (1979),  Janicka  et  al.  (1979),  Pope  (1982)].  The  two  partic¬ 
ipating  particles  in  the  above  example  will  have  the  new  values 

=  (1  -  «)?1  +  |a(§,  +  I2)  (15a) 

=  (1  -  a)92  +  |a(li  +  I2)  (15b) 

where  the  weight  a  is  varied  randomly  between  zero  and  one  according  to  the  flat  pdf  P(a)  =  1:  note  that 
other  choices  for  P(a)  could  have  been  made.  With  a  =0  the  particles  do  not  mix,  and  with  a=1  the  Curl 
model  is  reproduced.  To  produce  the  correct  decay  rate,  Nmix  must  be 

N™  =  3  Np  0.)  At  (15c) 

Once  again,  the  post- mixing  mass  fractions  and  temperature  of  each  particle  in  the  ensemble  are  ad¬ 
vanced  in  time  by  integration  of  the  chemical  kinetics  equations. 

The  modified  Curl  model  addresses  some  of  the  weaknesses  of  the  original  Curl  model;  however. 
Pope  (1982)  demonstrated  that  although  it  does  indeed  yield  a  continuous  pdf  when  starting  with  a  distribu¬ 
tion  which  is  composed  of  delta  functions  (provided  that  P(a)  is  continuous),  that  pdf  is  not  the  expected 
Gaussian.  The  peak  and  the  tails  of  the  pdf  exceed  those  of  the  Gaussian. 


(iii)  The  IBM  model  relates  each  particle  to  the  ensemble,  rather  than  to  a  partner.  The  equations  were 
given  above  and  will  not  be  repeated  here.  The  IBM  model  has  a  characteristic  which  is  evident  from  its 
parent  pdf  evolution  equation.  For  a  conserved  scalar  in  a  homogeneous  system  without  inflow  or  outflow, 
the  pdf  P(?,t)  evolves  according  to 


where  |  is  the  mean;  the  mean  does 


f -^l(M)Pl  =  0 

not  change  in  this  situation.  The  general  solution  is 


(16a) 


P(i,t)  =  +  ln_(^^)] 


(16b) 


where  fn  Is  a  function  defined  by  the  initial  condition  P(^,0).  It  is  clear  from  Eq.  16  (b)  that  the  lEM  model 
relaxes  the  pdf  (of  a  conserved  scalar)  in  a  shape-preserving  manner  and  does  not  approach  the  expected 
Gaussian  in  the  general  case. 


All  three  mixing  models  fall  to  yield  the  Gaussian  expected  of  a  conserved  scalar  pdf  in  the  limit  of 
large  time;  however,  as  mentioned  above,  the  significance  of  this  failure  is  not  obvious  in  the  reacting  case 
and  motivates  the  present  work. 

First  we  confirm  that  the  mixing  models  agree  identically  on  the  decay  of  the  variance  of  a  conserved 
scalar,  in  a  system  without  Inflow  or  outflow.  Consider  the  evolution  of  the  pdf  of  a  conserved  scalar,  P(^,t). 
Let  the  mixing  occur  at  a  rate  such  that  the  variance  decays  per  da^/dt=— 20)0^  (i.e.,  the  standard  devi¬ 
ation  o  decays  according  to  da/dt=  -(oo).  Each  of  the  three  ensembles  is  marched  forward  in  time,  without 
chemical  reactions. 
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Initially  symmetric  and  asymmetric  pdf’s  were  examined.  Figure  17  shows  typical  results  for  the  pre¬ 
scribed  decay  rate  a)=1000  s“'*.  There  are  no  discernible  differences  between  the  models,  all  of  which  re¬ 
produce  the  1000  s“^  decay  rate.  We  conclude  that  the  three  mixing  models  are  equivalent  (up  to  the  vari¬ 
ance)  for  the  case  of  a  conserved  scalar.  This  calculation  is  merely  a  test  of  the  present  implementation; 
the  results  of  Fig.  17  are  expected  [Pope  (1982)]. 

Next,  inflow,  outflow  and  chemical  reactions  are  reactivated  in  the  PaSR. 

The  fuel  is  taken  to  be  a  mixture  of  50%CO/50%H2  (by  vol .) .  The  kinetic  scheme  consists  of  1 1  species 
and  23  reactions  (Table  III).  The  inlet  conditions  are  1  atm  and  900K.  The  stoichiometry  of  the  premixed 
inflow  leads  to  a  PSR  temperature  of  1 740K  at  a  residence  time  of  5  ms  (1  BOOK  in  equilibrium) ,  but  to  blowout 
in  a  PFR.  The  PaSR  mixing  frequency  was  varied  in  the  range  10  -  10“^  Hz,  by  factors  of  it  has  been 
shown  that  this  range  more  than  covers  the  mixing  frequencies  encountered  in  a  practical  combustor  [Cor¬ 
rea  and  Braaten  (1993)].  The  computations  are  run  keeping  the  PaSR  volume  and  the  mass  flow  rate  fixed, 
so  the  residence  time  t  varies  with  the  mean  density  according  to  t  =  gV/m,  where  g  is  the  mean  density, 
V  is  the  reactor  volume  and  m  is  the  mass  flow  rate. 

In  the  interest  of  brevity,  only  the  temperature  will  be  discussed  here.  Figure  18  shows  the  evolution 
of  the  mean  temperature.  Results  are  shown  for  two  mixing  frequencies,  which  were  chosen  to  differentiate 
between  the  behavior  of  the  models  as  blowout  is  approached.  The  mean  temperature  attains  a  stochasti¬ 
cally  steady  state.  It  is  clear  that  the  IBM  model  sustains  combustion  at  both  the  chosen  frequencies,  where¬ 
as  the  Curl  and  modified  Curl  models  approach  blowout  at  the  lower  frequency.  The  initial  evolution  is  very 
similar  for  all  three  models:  the  PaSR,  which  was  initialized  with  particles  whose  composition  and  tempera¬ 
ture  were  set  to  the  PSR  results,  tends  to  cool  down  as  the  initial  particles  are  replaced  by  incoming  par¬ 
ticles.  The  Curl  and  modified  Curl  models  exhibit  a  greater  degree  of  fluctuations,  which  is  not  surprising 
as  blowout  is  approached. 

It  should  be  noted  that  we  present  all  frequencies  including  those  which  are  lower  than  the  s500 
Hz  expected  of  practical  turbulent  combustion.  The  models  become  increasingly  similar  in  their  conver¬ 
gence  paths  at  higher  frequencies,  although  the  Curl  and  modified  Curl  models  continue  to  exhibit  a  greater 
degree  of  fluctuations  in  the  mean. 

The  pdf’s  of  particle  age,  computed  from  the  steady-state  ensembles  that  were  obtained  using  each 
of  the  three  models  at  the  two  different  frequencies  which  will  be  examined  in  detail  below,  are  shown  in 
the  form  of  histograms  in  Fig.  19.  Also  shown  is  the  theoretical  pdf, 

P(i|.)  =1  exp(-i[./T)  (17) 

where  the  variable  is  the  age.  All  three  models  reproduce  the  theoretical  pdf  fairly  well;  the  deviations  are 
caused  by  the  finite  number  of  particles  in  the  ensemble  (=500). 


Table  III.  Kinetic  Mechanism  Used  for  C0/H2-Air  Combustion. 
(Forward  rate  constant  kf  =  A  moles-cm-s-K;  E  units:  cal/mole) 


No.  Reaction 

A 

b 

E 

.  1.  C0+0+M=C02+M 

3.20E+13 

0.0 

-4200.0 

2.  C0+0H=C02+H 

1.51E+07 

1.3 

-758.0 

3.  C04-02=C02+0 

1.60E+13 

0.0 

41000.0 

4.  H02+C0=C02+0H 

5.80E4-13 

0.0 

22934.0 

5.  H2+02=20H 

1.70E+13 

0.0 

47780.0 

6.  0H  +  H2=H20+H 

1.17E+09 

1.3 

3626.0 

7.  H+02=0H  +  0 

2.00E+14 

0.0 

16800.0 

8.  0+H2=0H  +  H 

1.80E+10 

1.0 

8826.0 

9.  H+02  +  M  =  H02  +  M 

2.10E+18 

-1.0 

0.0 

H20  enhanced  by  21 .0, 

C02  enhanced  by  5.0, 

H2  enhanced  by  3.3, 

CO  enhanced  by  2.0,  02  enhanced  by  0.0,  N2  enhanced  by  0.0 

10.  H  +  02  +  02  =  H02  +  02 

6.70E+19 

-1.4 

0.0 

11.  H+02  +  N2  =  H02  +  N2 

6.70E  +  19 

-1.4 

0.0 

12.  0H  +  H02  =  H20+02 

5.00E+13 

0.0 

1000.0 

13.  H+H02=20H 

2.50E  +  14 

0.0 

1900.0 

14.  0  +  H02=02+0H 

4.80E  +  13 

0.0 

1000.0 

15.  20H=0  +  H20 

6.00E+08 

1.3 

0.0 

16.  H2+M=H+H  +  M 

2.23E+12 

0.5 

92600.0 

H20  Enhanced  by  6.0.  H  enhanced  by  2.0,  H2  enhanced  by  3.0 

17.  02+M=0+0  +  M 

1.85E+11 

0.5 

95560.0 

18.  H+0H+M  =  H20+M 

7.50E+23 

-2.6 

0.0 

H20  enhanced  by  21 .0 

19.  H+H02=H2+02 

2.50E  +  13 

0.0 

700.0 

20.  H02+H02=H202+02 

2.00E  +  12 

0.0 

0.0 

21.  H202+M=0H  +  0H  +  M 

1.30E+17 

0.0 

45500.0 

22.  H202  +  H  =  H02  +  H2 

1.60E  +  12 

0.0 

3800.0 

23.  H202+0H  =  H20  +  H02 

1.00E  +  13 

0.0 

1800.0 

The  mean  temperature  is  displayed  against  the  mixing  frequency  in  Fig.  20.  At  high  enough  frequen¬ 
cies  (>1000  Hz),  all  three  models  approach  the  PSR  values.  At  lower  frequencies,  however,  the  pair-ex¬ 
change  models  tend  to  blow  out  while  the  IBM  model  sustains  combustion.  On  the  basis  of  the  results  for 
the  means,  it  may  be  provisionally  concluded  that  the  models  are  essentially  similar  at  the  high  mixing  fre¬ 
quencies  of  practical  interest.  This  conclusion  were  tested  in  greater  detail  by  examining  pdf’s  and  scatter- 
plots. 

The  temperature  pdf,  obtained  from  the  steady-state  ensembles,  shows  a  greater  degree  of  bi¬ 
modality  between  the  inlet  and  near-PSR  temperatures  with  the  pair-exchange  models  than  with  the  IBM 
model  (Fig.  21 ) .  These  differences  increase  as  the  frequency  drops  below  1 00  Hz.  New  particles  are  clearly 
more  quenched  by  the  former  models’  sequential  treatment  of  mixing  and  reaction  than  by  the  simulta¬ 
neous  treatment  in  the  IBM  model,  accounting  for  the  lower  mean  temperatures  encountered  above.  We 
also  note  that  the  equations  for  the  IBM  model  resemble  those  of  the  stirred  reactor,  which  admits  abrupt 
transitions  (bifurcations)  between  lit  and  unlit  states  in  high  activation -energy  systems.  Hence  it  is  not  sur¬ 
prising  that  the  IBM  model  approaches  blowout  in  a  more  abrupt  manner  than  the  pair-exchange  models. 
The  three  mixing  models  yield  similar  pdfs  at  the  high  temperature  end. 

These  results  indicate  that  the  choice  of  mixing  model  is  not  critical  in  the  distributed  reaction  regime 
of  lean  premixed  combustion,  as  long  as  the  turbulent  mixing  frequencies  are  above  about  1 000  Hz.  Three- 
dimensional  CFD  calculations  of  practical  high -intensity  combustors  indicate  in -flame  mixing  frequencies 
of  1000  Hz  and  above.  Hence,  in  such  combustors,  the  differences  between  models  are  ovenwhelmed  and 
blurred  by  the  finite  rates  of  the  chemical  reactions.  The  computational  advantages  of  the  IBM  model,  which 
is  well  suited  to  parallel  computer  architectures  whether  in  the  form  of  standalone  machines  or  a  networked 
cluster  of  workstations  operating  as  a  “virtual”  parallel  machine,  make  it  the  model  of  choice. 


TIMB,  s. 

Figure  17.  The  decay  of  the  variance  of  a  conserved  scalar,  computed  using  the  three  mixing  models. 

500  particles,  a)=1000  Hz,  At=0.1/co. 
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Figure  20.  Variation  of  mean  temperature  with  mixing  frequency. 

The  PaSR  has  shown  that  as  blowout  is  approached,  differences  arise  between  the  above  mixing 
models  which  are  commonly  used  in  CFD.  Which  is  “right”?  First,  are  PaSR  results  even  of  any  general 
significance? 

To  put  the  PaSR  in  context,  it  is  worthwhile  to  repeat  one  possible  vision  of  a  future  combustion 
model.  This  future  model  will  employ  the  Monte  Carlo  particle  tracking/pdf  transport  approach  to  turbu¬ 
lence-chemistry  interactions  and  be  dynamically  iterated  with  a  complex-geometry  mean  flow/two- 
equation  turbulence  CFD  model,  as  in  Section  2.1.  However,  it  will  not  use  reduced  chemistry  but  instead 
a  full  scheme  with  a  large  number  of  species  (eventually,  as  many  as  50).  This  degree  of  precision  in  the 
chemistry  is  required  to  address  the  issues  of  the  day  (NO^  down  to  the  level  of  5  ppm,  CO,  unburned  hydro¬ 
carbons,  flame  stability  under  extreme  conditions,  and  so  on).  Within  each  computational  cell  the  system 
is  spatially  homogeneous  but  unmixed,  as  in  the  present  PaSR.  Since  full  chemistry  is  too  expensive  to  use 
in  calculations  of  general  inhomogeneous  flows  at  present,  here  we  have  in  effect  concentrated  on  turbu¬ 
lence-chemistry  interactions  in  the  single  cell  (i,e,,  on  a  single  PaSR).  Hence,  data  from  inhomogeneous 
systems  are  not  directly  usable.  Experimental  data  on  a  PaSR  are  needed  to  decide  which  model  is  most 
"right”  physically  and  to  make  further  choices  on  mixing  closure. 
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ure  21 .  Pdf  of  temperature  at  three  selected  mixing  frequencies. 


Section  3 
CONCLUSIONS 

We  have  conducted  experiments  on  bluff- body  stabilized  flames  of  CO/H2  and  CH4/air.  By  correct¬ 
ing  for  fluorescence,  we  have  made  Raman  measurements  in  bluff-body  CH4-air  flames;  however,  the 
errors  in  certain  species  (e.g.,  CO  and  CO2)  may  be  so  large  that  models  should  not  be  changed  on  the 
basis  of  such  data  alone.  This  of  course  detracts  severely  from  the  utility  of  expensive  laser-based  mea¬ 
surement  techniques.  In  refereed  publications  and  in  oral  presentations,  we  have  raised  the  consciousness 
of  the  university  and  industrial  research  community  to  this  issue. 

We  have  examined  modeling  approaches  of  varying  complexity. 

1 .  The  assumed  shape  pdf  model:  it  was  dismissed  because  it  is  not  extendable  to  multi-scalar  chemistry. 

2.  The  “conditional  moment  closure”  model:  it  was  dismissed  because  it  is  applicable  only  when  the 
fluctuations  about  the  conditional  mean  are  small  (i.e.,  as  in  1^-intensity  turbulent  combustion). 

3.  The  scalar  pdf  model,  which  does  not  treat  the  velocity  part  of  the  pdf.  Instead  this  model  uses  standard 
turbulence  modeling  to  supply  the  scalar  (and  momentum)  fluxes.  It  is  a  reasonable  engineering  tool  in 
3D  design,  but  it  is  not  state-of-the-art  in  terms  of  research  and  hence  cannot  be  used  to  develop 
advanced  tools  of  the  future. 

4.  The  joint  velocity-composition  pdf  transport  model.  During  this  contract  period,  we  have  extended  joint 
pdf  modeling  from  parabolic  (jet)  flames  to  elliptic  (bluff-  body  recirculation-stabilized)  flames  as  found 
in  practical  burners.  The  axisymmetric  elliptic  mean  flow  CFD  model  and  the  pdf  model  communicate 
with  each  other  iteratively.  The  speedup  achieved  in  parallelizing  particle  tracking  pdf  computations  - 
accomplished  in  the  present  reporting  period  of  this  contract  -  shows  that  the  pdf  computations  are 
tractable  in  practical  design  codes. 

The  Monte  Carlo  joint  PDF/CFD  model  has  been  assessed  by  comparisons  with  Raman  data  on 
means  and  rms  of  mixture  fraction,  major  species  and  temperature  in  non-premixed  bluff-body  stabilized 
CO/H2  and  CH4  flames.  For  brevity  and  because  most  of  the  points  were  made  in  the  subsequent  CH4 
flame,  the  CO/H2  flames  were  not  discussed  here.  Details  are  available  in  Correa  and  Pope  (1992).  The  CH4 
flame  was  discussed  in  detail,  demonstrating  several  points: 

1 .  The  bluff  body  burner  provides  a  strongly  turbulent  field  leading  to  localized  extinction,  without  the  need 
for  a  pilot  flame.  Thus  the  two-stream  nature  of  the  problem  is  preserved,  unlike  many  piloted  jet  flame 
studies  where  the  composition  or  the  excess  enthalpy  of  the  pilot  flame  can  cause  modeling  difficulties. 
The  present  recirculation-stabilized  flame  is  also  much  closer  to  practical  burners. 

2.  Given  the  similarity  in  scatterplots,  it  is  clear  that  the  pointwise  structure  of  the  above  bluff  body  flame  and 
the  piloted  jet  flame  are  quite  similar.  A  greater  degree  of  local  extinction  is  measured  here. 

3.  The  limitations  of  pdf  shape  assumption,  statistical  independence  of  scalars,  and  gradient  diffusion,  are 
removed  in  this  joint  pdf  model.  The  consequences  can  be  seen  in  joint  scatterplots  and  in  convective 
fluxes. 

4.  The  calculated  CO  scatterplots  have  maxima  of  3%  whereas  the  data  peak  at  about  1 0%,  well  above  the 
flamelet  maxima.  Similar  10%  levels  were  measured  another  bluff  body  flame  and  In  pilot-stabilized 
flames,  and  2-3%  peaks  were  predicted  in  the  latter  using  the  4-step  scheme.  Hence  this  discrepancy 
on  CO  maxima  has  appeared  in  diverse  circumstances  (but  always  in  combustion  gases  that  are  near 
local  extinction).  There  are  many  potential  contributors,  including:  (i)  the  assumption  of  a  steady  state  for 
the  radicals  in  the  4-step  mechanism;  (ii)  the  errors  in  Raman-based  CO  and  CO2  data;  and  (iii)  neglect 
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of  phenomena  such  as  unsteady  flamelets  or  micromixed  gases  (perfectly  stirred  reactors)  which  have 
been  shown  to  lead  to  high  CO.  In  intense  turbulence,  however,  the  microscale  may  be  better  simulated 
by  spatially  homogeneous  systems. 

It  seems  clear  that  CO  is  an  important  clue  to  the  microstructure  of  highly  turbulent  combusting 
gases,  and  that  accurate  measurements  will  be  critical. 

The  acquisition  of  Raman  data  and  the  3-velocity/5-scalar  joint  pdf  calculation  in  this  bluff  body 
methane  flame  took  each  “discipline”  to  the  limits  of  the  present  state  of  the  art.  Any  significant  further  prog¬ 
ress  is  likely  to  require  improvements  in  major  species  measurements,  complementary  velocity  and  minor 
species  measurements,  and  parallel  computers.  Reduced  chemistry  schemes  which  relax  steady-state 
assumptions  (and  are  likely  to  require  additional  scalars)  will  have  to  be  developed  and  assessed  in  simpler 
contexts. 

We  have  shown  that  an  intense  combustion  process  should  not  be  viewed  as  the  motion  of  a  high¬ 
speed  flamefront;  as  the  turbulence  amplitude  and  frequency  increase,  the  relevance  of  reaction-diffusion 
interfaces  (flamefronts)  decreases.  Ultimately,  the  velocity  field  and  the  spatial  structure  become  less  impor¬ 
tant  and  a  stirred  reactor  description  would  be  more  appropriate  at  the  micro-scale  level. 

To  the  latter  end,  we  developed  the  “PaSR”  or  partially  stirred  reactor  model  for  spatially  homoge¬ 
neous  (but  not  mixed)  reacting  flows.  The  PaSR  offers  an  alternative  to  traditional  evaluations  in  laminar 
flames,  which  are  not  relevant  to  turbulent  combustion. 

We  have  demonstrated  that  particle  tracking  pdf  calculations  are  easily  parallelized,  whereas  flow 
codes  “prefer”  vector  computers.  The  parallel  implementation  was  targeted  at  distributed  memory  MIMD 
(multiple  instruction,  multiple  data  set)  computers  such  as  the  Intel  iPSC/860  and  Intel  Delta.  From  the  huge 
speed-up  (factors  of  100)  in  execution  time,  it  is  clear  that  such  computations  are  indeed  practical  in  the 
networked— workstation -intensive  environment  of  the  typical  gas-turbine  design  community.  It  was  also 
clear  that  distributed  memory  MIMD  parallel  architectures  are  well  suited  to  the  particle -tracking  PaSR  al¬ 
gorithm.  This  remains  true  so  long  as  the  turbulence  model  does  not  require  continuous  communication 
between  “particles.”  Thus,  the  lEM  model  is  preferable  to  pair-exchange  models  on  parallel  or  networked 
computers. 

We  have  used  the  PaSR  model  in  first-ever  studies  of  a  full  (27-species/77- reactions)  methane 
oxidation  scheme  interacting  with  prescribed  turbulence.  The  calculation  involved  600  to  700  particles  and 
therefore  —20,000  o.d.e.’s.  The  PaSR  smoothly  merges  into  PSR  and  PFR  limits  as  mixing  frequency  be¬ 
comes  large  or  small,  respectively.  The  degree  of  turbulence  in  the  time-  histories  of  species  is  apparent, 
e.g.,  greatest  for  CO  (among  major  species)  as  blowout  is  approached.  The  “skeletal”  mechanism  used 
as  the  basis  for  the  4-step  reduced  scheme  in  laminar  flames  cannot  reproduce  the  full  scheme’s  results 
at  low  frequency,  where  ignition  chemistry  becomes  important. 

We  have  also  applied  the  PaSR  to  assess  a  3-variable  (in  the  non-premixed  case)  reduced  chemis¬ 
try  scheme  for  a  hypothetical  heavy  hydrocarbon.  The  agreement  between  the  starting  scheme  and  the 
reduced  scheme  was  excellent  on  the  mean  temperature  and  on  the  mean  fuel  mass  fraction.  Fuel  pyrolysis 
is  explicitly  recognized  as  a  degree-of-freedom  in  the  reduced  scheme.  The  latter  is  a  desirable  feature. 
The  agreement  on  mean  Y*,  which  is  the  combined  variable  based  on  the  CO  and  H2  mole  numbers,  was 
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not  as  good.  Potential  contributors  to  this  discrepancy  have  been  identified.  On  balance,  the  reduced 
scheme  exhibits  encouraging  performance. 

We  have  also  used  the  PaSR  to  directly  compare  the  lEM  model,  the  Curl  model,  and  the  modified 
Curl  model  in  the  context  of  a  full  scheme  in  prescribed  turbulence.  The  results  indicate  that  the  choice  of 
mixing  model  is  not  critical  in  the  distributed  reaction  regime  of  lean  premixed  combustion,  so  long  as  the 
turbulent  mixing  frequencies  are  above  about  1000  Hz.  Three-dimensional  CFD  calculations  of  practical 
gas-turbine  combustors  indicate  in-flame  mixing  frequencies  of  1 000  Hz  and  above.  Hence,  in  such  com¬ 
bustors,  the  differences  between  mixing  models  are  overwhelmed  and  blurred  by  the  finite  rates  of  the 
chemical  reactions.  The  computational  advantages  of  the  lEM  model,  which  is  well  suited  to  parallel  com¬ 
puter  architectures  whether  in  the  form  of  standalone  machines  or  a  networked  cluster  of  workstations  op¬ 
erating  as  a  “virtual”  parallel  machine,  make  it  the  model  of  choice. 

We  have  only  briefly  summarized  the  development  and  some  key  conclusions  of  the  PaSR  model. 
Many  more  results  and  sensitivity  analyses  are  given  in  our  refereed  literature  (Section  5).  The  PaSR  model 
has  also  been  adopted  as  an  investigative  tool  by  other  researchers,  e.g.,  at  the  University  of  Sydney  and 
at  UC-Berkeley. 
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Section  5 

WRITTEN  PUBLICATIONS 


This  report  is  not  comprehensive  in  terms  of  the  details  of  the  work  that  was  performed.  Additional 
description  may  be  found  in  the  several  refereed  reports  that  were  produced  as  a  result  of  this  contract. 

Refereed  publications  based  on  work  performed  in  and  written  during  the  reporting  period  were: 

1 .  Correa,  S.M.  and  Pope,  S.B., “Comparison  of  a  Monte-Carlo  PDF/Finite-Volume  Mean  Flow  Model  with 
Bluff-Body  Raman  Data,”  Twenty-Fourth  Symposium  (International)  on  Combustion,  The  Combustion 
Institute,  Pittsburgh,  PA,  pp.  279-285,  1992. 

2.  Correa,  S.M. .“Turbulence-Chemistry  Interactions  in  the  Intermediate  Regime  of  Premixed  Combus¬ 
tion,”  Comb,  and  Flame,  93,  pp.  41  -60,  1993. 

3.  Correa,  S.M.,  and  Braaten,  M.E., “Parallel  Simulations  of  Partially- Stirred  Methane  Combustion,”  GE 
Class  I  Report  92CRD273,  November  1992;  Comb,  and  Flame,  94,  pp,  469-486,  1993. 

4.  Correa,  S.M., “Models  for  High-Intensity  Turbulent  Combustion,”  Computing  Systems  in  Engineering, 
Vol.  5,  No.  2,  pp,  135-145,  1994. 

5.  Correa,  S.M. .“Assessment  of  a  3— Variable  Reduced  Kinetic  Scheme  in  Prescribed  Turbulence,”  in  press 
J.  Prop.  Power,  1994. 

6.  Correa,  S.M.,  Gulati,  A.,  and  Pope,  S.B.. “Raman  Measurements  and  Joint  PDF  Modeling  of  a  Non-Pre- 
mixed  Bluff  Body— Stabilized  Methane  Flame,”  Twenty-Fifth  Symposium  (International)  on  Combustion, 
Irvine,  CA,  July  31  -  August  5,  1994. 

7.  Correa,  S.M.  and  Dean,  A. J., “Turbulent  Broadening  of  Autoignition  Limits,”  Twenty— Fifth  Symposium 
(International)  on  Combustion,  Irvine,  CA,  July  31  -  August  5,  1994. 

8.  Correa,  S.M.,“A  Direct  Comparison  of  Pair-Exchange  and  lEM  Models  in  Premixed  Combustion,”  sub¬ 
mitted  to  Comb,  and  Flame,  April  1994. 


Section  6 

PROFESSIONAL  PERSONNEL 


Dr.  Sanjay  M.  Correa  (PI),  Dr.  Mark  E.  Braaten,  and  Ms.  Janet  Sober  worked  on  this  program.  Dr.  Iris 
Hu  became  a  member  of  the  team  in  the  last  few  months  of  the  program. 


Section  7 

INTERACTIONS/TECHNOLOGY  TRANSFER 

The  computational  technology  developed  in  the  program  has  been  widely  disseminated  in  the  archi¬ 
val  literature  and  in  invited  talks  and  seminars,  and  elements  of  it  are  used  in  3D  CFD  models  of  combustors. 
The  measurements  provide  corroborating  data  not  only  for  the  models  developed  here  but  also  for  others. 
The  research  has  also  provided  techniques  for  laser  probing  of  practical  combustors,  e.g.,  multi-cup  sec¬ 
tors  are  now  being  probed  with  Raman  spectroscopy.  Since  we  have  shown  that  the  computational  meth¬ 
ods  developed  in  this  research  program  can  be  fielded  within  the  CFD  codes  used  in  design,  ongoing  trans¬ 
fer  of  the  new  research  is  assured. 


Section  8 
INVENTIONS 

There  were  no  inventions  in  the  reporting  period. 
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